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Foreword 


The  first  seven  sections  of  this  report  contain  brief,  but 
technical  descriptions  of  research  on  numerical  fluid  mechanics 
that  has  been  carried  out  at  Brown  University  under  contract 
with  the  Army  Research  Office.  More  detailed  descriptions  of 
parts  of  that  research  are  provided  in  six  appendices. 

The  sections,  except  Section  8,  are  arranged  in  chronological 
order  according  to  when  we  began  work  on  the  various  problems. 

-  At  the  start,  three  and  a  half  years  ago,  the  finite  fluid 
element  method  (Section  1)  was  the  only  one  under  consideration. 
Things  went  badly  with  that  method,  and  progress  with  its 
implementation  went  far  more  slowly  than  we  ever  anticipated. 

As  a  result,  by  the  end  of  the  first  year,  a  second  method,  which 
originally  was  developed  for  a  check  on  results  of  the  first, 
had  become  by  far  the  more  promising  of  the  two.  Most  of  this 
report  (Sections  3-6)  describes  progress  we  have  made  with  the 
application  of  biased  differences  (Section  2)  to  a  variety  of 
fairly  difficult  problems  of  numerical  fluid  mechanics.  Finally, 
in  the  last  six  months  of  the  period  covered  by  this  report, 
the  major  difficulties  with  finite  fluid  elements  were  overcome, 
so  it  became  possible  to  begin  a  comparison  of  the  two  methods 
(Section  6).  Preliminary  indications  are  that  both  methods  are 
reliable,  and  both  are  considerably  more  efficient  than  a  third 
method  with  which  they  have  been  compared. 

Section  7  is  a  report  of  progress  with  a  boundary  integral 
method  that  is  not  closely  related  to  the  others  except  by  being 
numerical,  and  Section  8  is  a  description  of  the  kinds  of  graphical 
software  we  had  to  develop  for  interpretation  of  our  numerical 
computations . 

The  research  reported  here  has  been  carried  out  by  F.  E.  Bisshopp, 
R.  B.  Caswell  (principal  investigators),  E.  W.  Fieri,  M.  E.  Michaud, 
and  T.  G.  McKee  (research  assistants  in  Applied  Mathematics). 
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1.  Finite  fluid  elements. 

As  initially  formulated,  the  finite  fluid  element  method  dealt 
with  localized  distributions  of  mass,  momentum  and  energy.  The 
mass  density  and  momentum  density  of  a  fluid  were  represented  as 

N 

p(x,t)  -  I  mi(t)f1(x  -  Xjl(  t ) ) 

N 

pu(x,t)  ~  £  mix1f iC*  ~  xi) 


with 


distribution  functions  in  n 

|x-X 


fi  " 


2/  2 
/al 


=  1, 


?  or  3  dimensions. 


The  equations  of  fluid  mechanics  were  then  used,  along  with  least 
square  fitting  of  the  approximations,  to  obtain  ordinary  differential 
equations  for  the  fluid  element  parameters  m1(t),  X^t)  and  a^(t). 

A  detailed  description  of  the  method  is  included  in  appendix  A: 
here  we  will  report  on  two  major  difficulties  and  what  has  been 
developed  to  overcome  them. 

In  the  course  of  our  attempts  to  implement  a  finite  fluid  element 
method  we  found  first  that  the  original  formulation  was  simply  too 
complicated  to  be  of  practical  value  —  it  was  never  successfully 
employed  for  a  two-dimensional  flow.  Instead,  two-dimensional 
flow  has  been  treated  by  an  algorithm  that  is  based  upon  the 
original  one,  but  is  greatly  simplified  as  follows: 

Simulation  of  the  equation  of  continuity  has  been  dropped  in 
favor  of  fixed  particle  masses,  thus  bringing  the  method 
somewhat  nearer  to  a  particle-in-cell  method.  The  basic 
approximation  is 

N 

P ( x, t )  ~  £  nijfj/x  -  X1(t)) 

where  the  integral  of  ft  over  all  of  1,  2  or  3-dimensional  space  is 


<ft>  =  I- 


W&’ 


2 

L  U 

The  mass  density  at  the  center  of  the  1°  element  Is  assigned 
the  (approximate)  value 


N 


mj<ftf ^ 


P 


i 


where 


<fiV 


<’(°i 


°j))n/2 


and  the  parameters  a.,  that  characterize  element  diameters  are 
adjusted  to  give  an  overlap  of  the  i  element  with  2,  6  or  12 
neighbors  in  1,  2  or  3  dimensions.  Some  Indication  of  ways  to 
assign  the  a's  is  given  in  appendix  A,  and  further  investigation 
of  the  effect  of  that  choice  on  accuracy  of  the  method  is  still 
in  progress. 

Motion  of  the  centers  of  the  elements  is  governed  by  the 
approximation  of  the  fluid  momentum  equation. 


Pi *  =  P1£(£i)-(.2P)i  +  (V‘t)1 

for  viscous  flow,  or 


PiX’i  =  Pl£(Xi)  -  (vp^  - 

for  flow  in  a  porous  medium.  The  body  force  £(x)  has  presented 
no  difficulty,  and,  for  a  barotropic  fluid,  it  has  been  found 
that  an  adequate  approximation  of  the  pressure  gradient  is 

(_ZP)i  »  P'(  Pl)(  V  p)j_ 

Orp>i  -  j  V<fiV 

N  m . 

=  ?  I  -??  (X1  -  X.)  • 

1  c,+o,  '  J 


where 
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A  second  difficulty  with  the  initial  formulation  was  found  In 
the  attempt  to  approximate  the  effect  of  viscous  stress  The 
quadratic  fitting  of  velocities  outlined  in  appendix  A  gave  a 
rather  poor  estimate  of  the  second  gradients  that  appear  in  a 
direct  evaluation  of  (V-t)^.  It  has  been  found  now  that  a 
much  better  representation  of  viscous  effects  can  be  obtained 
in  terms  of  two  estimates  of  first  gradients.  In  the  simplified 
algorithm  for  viscous  flow  there  is  a  viscous  stress  field 

N 

T(x,t)  ~  l  t1fi(x  -  X1(t) )  , 


and  its  divergence  is  approximated,  similarly  to  _Vp  ,  as 

N  (XrX-)-t- 

(jz-  r),  -  ?  i  i  s3  «,iiy  • 

1  1  of  +  07  J 


The  coefficients  t^  are  defined  implicitly  by  approximate  values 
of  the  stress  field  as 

1  <fiV 


where 

=  x(  p4)(  v-u)1i  +  y(  Pj_)(  (_v  u)^  +  (_v  u)^) 

and  ( V  u)^  is  the  outer  product 

N  {X.-XJ  . 

(-  M)1  =  2  l  0^+  a?  -J  <fif P  ' 

1  °i  j 

Investigation  of  performance  of  the  1  and  2-dimensional  versions 
of  the  simplified  algorithm  is  still  in  progress.  Results  on  sound 
waves  and  shock  waves  will  appear  in  the  Ph.D.  dissertaion  of 
E.  W.  Fieri,  and  he  is  scheduled  to  present  a  poster  session  on 
the  subject  at  the  SIAM  Meeting,  Alexandria,  Va . ,  5-7  June  19^0 
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2.  Biased  differences. 

The  method  of  biased  differences  was  originally  formulated  and 
developed  as  a  standby  --  it  was  to  be  used  for  independent 
checks  of  results  of  the  finite  fluid  element  method.  As  it 
turned  out,  however,  very  few  problems  with  it  were  encountered, 
and  its  present  state  of  development  is  considerably  advanced 
beyond  that  of  the  finite  fluid  element  method.  It  has  been 
applied,  with  promising  results,  to  relatively  difficult  problems 
of  2  and  3-dimensional  unsteady  flow.  In  this  section,  we  give 
a  detailed  description  of  the  application  of  the  method  to  the 
simple  problem  of  Burgers'  equation  in  one  dimension 
Given 

ut  +  uux  -  vuxx" 

the  left-hand-side  is  a  material  derivative,  i.e.  the  total 
time  derivative  along  a  particle  path  defined  by  X  =  u(X(t),t), 
and  the  mixed  Eulerian-Lagrangian  formulation  of  the  Burgers' 
equation  is 

X  =  u 


By  contrast  with  finite  fluid  elements  that  move  along  particle 
paths,  biased  differences  are  defined  on  a  grid  that  is  fixed 
in  space  and  has  the  local  skeleton. 


There  are  several  possibilities  for  defining  the  bias  line 
(approximate  particle  path  in  this  case)  and  performing  the 
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time  integration  to  various  orders  of  accuracy.  By  experiment 
it  has  been  found  that  the  following  relatively  simple,  low  order 
scheme  gives  quite  reliable  results. 

1)  u  is  taken  to  be  the  linear  interpolation  of  vJor  v  ) 
c  __  X/ v  r 

and  vc  at  a  point  xc 


u. 


=  v 


vi 


(vn  >  0) 


u. 


=  V 


(x/r*c) 

h 


+  v 


( x 

'  c  c 


(vc  £  0) 


2)  The  bias  line  is  defined  by  the  forward  difference 

x  -  x  +  ku  ,  and  elimination  of  x  -x  gives 
c  c  c  c  c 


uc  - 


<vc  i  °> 


uc  - 


1  +  ft  VV0) 


(vc  £  o) 


3)  The  equation  for  u  is  approximated  by  a  backward  difference 
in  time  and  a  central  difference  in  space  to  give 


-  ?U  )  =  u 
r  c  c 


The  low  order  integration  scheme  has  a  local  truncation  error 

o 

of  0(k  )  in  its  time  integrations,  but  it  has  not  been  found 
necessary  to  improve  that  by  introduction  of  central  differences. 
An  indication  of  how  the  combined  forward  and  backward  differences 
behave  like  central  differences  is  given  in  appendix  B. 


o 


3.  Two-dimensional  unsteady  flow 

Introduction  of  body  forces  and  application  of  biased  differences 
to  2  and  3-dimensional  versions  of  Burgers'  equation  is  straight¬ 
forward  and  will  not  be  discussed  here.  We  turn  now  to  the  more 
difficult  problems  of  incompressible  fluid  mechanics.  In  mixed 
Eulerian-Lagrangian  form  the  equations  of  motion  are 

X  =  u 

u  =  f  -  Vir  +  vAu 

£‘U.  =  0 


where  f  is  a  specified  body  force,  tt  is  p/p  ,  and  v  is  the 
kinematic  viscosity. 

As  in  the  one-dimensional  example,  the  grid  is  fixed  in  space, 

u„  refers  to  a  central  node  at  time  t,  v  refers  to  the  same  node 
— c  — c 

at  time  t-k,  x  is  the  position  of  a  node,  and  x  is  the  position 
c  c 

at  time  t-k  of  the  bias  line  that  passes  through  xc  at  time  t. 


The  local  skeleton  (In  plan  views)  is 
— n 


— n 


H„ 


u 

— c 


HS 


He 


Hw 


( time  t  ) 


—s 


He 


(  time  t-k  ) 


The  forward  difference  approximation. 


x  =  x  +  ku„ 
c  c  — c 


and  linear  Interpolation  for  uQ  now  gives  simultaneous  linear 
equations  for  the  two  (or  three)  components  of  uc .  The  equations 
depend  upon  which  quadrant  (or  octant)  contains  xc  and  that  is 
decided  by  the  signs  of  the  components  of  v_c.  If  both  components 
of  v  are  positive,  for  example,  x  is  in  the  third  quadrant  and, 

v*  C 

with  subscripts  1  and  2  to  denote  components  of  u  and  v  , 


•  v  rv  WfrJvWsa  ~ 
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<1  +  £<vlc  -  vle>,ulc 


+  £<vlc  -  vl>?c  *  vlc 


ft' 


2c 


2e 


>“lc  +  <>  +  ft'SC 


23 


>>u2c  -  v?c 


This  gives  initial  values,  u,  for  the  integration  of  the  momentum 
equation. 

To  find  the  velocity  field,  the  momentum  equation,  is  first 
approximated  by  a  backward  time  difference  with  no  space  differences 
included  yet: 


Hc  =  +  k(f(xc)  -  (Vtt)c  +  v(Au)c)  . 

Then  the  central  difference  approximation  of  _V-u  =  0  =  AV/u 
gives  the  pressure  equation. 


'u 


+  tt„  +  ir  +  it..  -  4ir 


+ 


_  v,,fl(£e)  “  fl^-w^  .  f?^-n^  "  f?^s^ 

m  ?  +  p 

h  ^Ule  “  ulw  .  uln  "  ulsN 
k  '  2  2  ' 


With  values  of  -n  determined  by  approximate  solution  of  the  Poisson 
equation,  central  space  differences  in  the  time-integrated  momentum 
equation  finally  give 


ulc  '  h^uln  +  uls 

+  ule  + 

ulw 

- 

»«i«) 

— 

«lc 

krl(*c> 

,  k  *e 

-  IT 

W 

i 

+  h 

2 

u2c“  ^l^u2n  +  u2s 

+  u2e  + 

U2w 

- 

"u2c> 

U2c 

4- 

kf?(xc) 

k  ^n 

-  IT 

s 

T 

h 

2 

In  practice,  it  has  been  found  that  the  implicit  equations 
for  it  and  u  at  the  nodes  can  be  treated  quite  efficiently  by 
iterations,  with  starting  values  obtained  from  data  at  time  t-k. 

There  remains,  as  always,  the  treatment  of  boundary  conditions 
So  far,  only  rectangular  geometries  with  boundaries  parallel  to 
a  coordinate  axis  have  been  considered.  At  a  southern  boundary, 
for  example,  with  xc  on  the  boundary  V -u  =  0  implies 


u2(xic'  x2c  +  y)  ""  uo 


2 


and 


vAu 


2  ~  .  2  U2n  * 
h 


At  y  =  §-h  the  normal  component  of  the  momentum  equation  then 
gives  the  pressure  boundary  condition. 


*c  =  *u 


^(f^Xjj)  +  f2(xc))  h  u?  +  *  k^u?n 


u?n>' 


This  too  is  implicit,  and  iteration  has  proved  effective. 

Appendix  C  is  a  description  of  an  application  of  the  method 
to  flow  in  a  channel.  The  results  obtained  there  compared 
very  well  with  results  obtained  by  a  far  more  elaborate  finite 
element  method . 


4.  Thermal  convection 


At  this  time  the  only  three-dimensional  unsteady  flow  that  has 
been  simulated  is  thermal  convection.  The  Bousslnesq  approxima¬ 
tion  for  nearly  Incompressible  fluids  (i.e.  liquids)  is  governed 
by 

X  =  u 

u  =  -  _Vtt  -  a'T  -  <T>)£  +  vAu 
T  =  <AT  +  <J> 

V*u  =  0 

where 

ir  is  —  4-  V  - nearly  constant  density 

Po 

g  is  -  V V  --  gravitational  force  and  potential 
a  is  the  coefficient  of  thermal  expansion 
<T>  is  the  average  temperature 
v  is  the  kinematic  viscosity 
k  is  the  thermal  diffusivity 
the  viscous  dissipation,  4>  ,  is  usually  negligible. 

In  essence,  the  approximation  is  derived  by  neglecting 
compressibility  everywhere  but  in  the  part  of  buoyancy  that  is 
not  derivable  from  a  potential  . 

The  biased  difference  scheme  for  these  problems  is  a  forward 
difference  of  X  =  u  to  define  a  bias  line  followed  by  backward 
differences  for  the  momentum  and  energy  equations,  and  again 
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there  is  a  Poisson  equation  for  the  pressure.  Since  initial 
data  u  and  T  is  used,  it  is  preferable  to  solve  the  forward 
difference  equations  for  x  -  x  instead  of  u  .  With  directions 

V/  vj  _  C 

relative  to  central  nodes  defined  as 


a. 


y 


the  equations  for  displacements  when  the  components  of  _vc  are 
all  positive  are 

<X  +  K(vlc-vlw)>«l  *  K<vlc-vU>S2  +  H^vlc_vld^53  '  kvlc 


2c 


2w 


K<v3c-v3» 

where  6  is  x  -x  * 

—  — c  — c 

other  octants. 

Given  _6  as  above, 
is 


)6i  +  (1  +  h^v2c'v?s^6?  +  h^v?c'v?d^63  kV2c 


)61  +  k^v3c‘v3s^2  +  (1  +  h(v3c'v3d^63  ~  kv3c 


Similar  sets  of  equations  define  6  in  the 


the  initial  data  for  the  backward  differences 

(v0-vw)4  -  <ic-“s>Tr  -  f-VVir 
-  <VTB>Tr  -  <VTd$ 


and  the  backward  equations  are: 
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V  rf<Te  +  Tw  +  Tn  +  Ts  +  Tu  +  Td  -  6tc>  -  \ 


u  -  u  +  u  ,  +  u„  +  uD  +  u„  +  u .  -  6u„) 

— c  h2  — e  — w  ~n  — 8  — u  — a  — c 


=  u  +k[(Vir)  -  a£(T  -  <T>)  ] 


1  , 


^  (it  +  it  +ir  +  tt  +  it.  +tt,  -  6it  )  =  (V  *u)  -  a  (g-VT) 

2'e  w  n  s  u  d  c/v - 'c  v*  —  '< 


Some  results  on  thermal  Instability  of  a  fluid  heated  from 
below  (the  Benard  problem)  are  given  in  appendix  D. 
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5.  Two  phase  flow 

A  considerable  effort  was  devoted  to  the  formulation  of  a 
biased  difference  algorithm  for  flow  of  two  immiscible  liquids 
in  a  porous  medium.  In  its  full  generality  the  problem  there 
is  rather  complicated,  but  there  is  nothing  in  its  structure 
that  precludes  an  application  of  the  method  to  it.  For  the  most 
part,  however,  we  have  concentrated  on  a  simplified  version 
where  the  sum  of  the  volumetric  flow  rates  of  the  'wetting'  and 
'nonwetting'  phases 


is  specified  in  advance  to  be  an  irrotational ,  incompressible 
flow  field.  Then  the  ratio  of  the  proportions  of  void  volume 
filled  by  the  phases  (saturation) 

-  vw/(vn  +  Vw) 


is  governed  by  the  saturation  equation, 

||  +  F'(s)u-Vs  =  V*$(s)Vs. 

In  this  case  the  forward  difference  for  the  bias  line, 

*c  *  Ic  +  kP’(8(Ic^« c  ‘ 


is  not  so  easily  solved  as  before  because  the  derivative  of  the 
fractional  flow  rate,  F'(s),  is  not  a  linear  function  of  s.  In 
fact, 


c)  --  si1:') 


(s2+  (l-s)V 


1 8  more  typical  of  the  kind  of  problems  encountered  in  waterflooding 
of  oil  producing  reservoirs 
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The  case  where  F'(s)  is  given  as  above  and  /9  is  constant  has 
been  investigated  in  detail  by  Marion  C  Michaud  in  her  Ph  D 
dissertation.  The  dissertation  is  included  here  as  appendix 
E:  a  shorter  version  for  publication  is  being  prepared 


6.  Compressible  flow 


A  start  has  been  made  on  the  application  of  b'ased  differences 
to  compressible  flow.  In  the  simple  barotrooic  case 

P  =  PP 
X  =  u 

p  =  -  pV/u 

u  =  - ?Vp  +  -jj  u  +  ^  I*u  I). 

Comparisons  of  results  in  the  one  dimensional  case  have  been 
made  between  the  finite  fluid  element  method,  a  biased  central 
difference  algorithm,  and  the  biased  forward -backward  difference 
algorithm. 

The  central  difference  algorithm  is 


xc  +  ^k(pc  +  pc> 

-  (J-  -  gk(yu)c 

^  (1  +  i  £(ur-uf)) 

_  (  P„-  Po ) 

uc  =  luc  •  *k(<Vp)c  +  ~  gh  '  '  } 

+  jvk((Au)c  +  -^(u^+ur) )  ]/(l  +  -^4j) 


where  bars  denote  linear  interpolation  at  xc-  Iteration  is  carried 
out  on  all  three  equations  in  the  sequence  indicated  above. 

The  much  simpler  forward -backward  algorithm  is 


cc  ■  pe/(1  +  4(ur-ui>> 

uc°  ‘  'Pi)  +  %Vuf)1/(1 

n  “ 


in  which  only  the  last  two  equations  are  iterated. 


The  forward -back ward  finite  fluid  element  method  that  was 
used  ia 


m .  =1 


X4  = 


Xi  +  ku^ 


O  4  = 


Pi  - 


(vp)jl  = 


v  r(Xj  -  xi )2 

(x,-x.) 

2  V  /  T- 


°J+0i 


(Vu)1 


tj  <f1r1> 


(2  neighbors) 

(i  and  2  neighbors) 

(2  neighbors) 


„  (x.-x1) 

=  2  ['  g-1^  g —  Uj  <f1fj>  (2  neighbors) 


u,  = 


(vu)l  - 

5i ' 


t,  <f,f: 


°Jwi 


'J 


in  which  the  last  three  equations  are  iterated. 

In  simulations  of  one-dimensional  sound  and  shock  waves,  all 
three  methods  worked  well;  and  we  were  led  to  conclude  that 
central  time  differences  give  no  significant  improvements  over 
the  simpler  methods.  Computation  times  for  the  forward -backward 
biased  difference  and  finite  fluid  element  methods  were  comparable 
and  significantly  shorter  than  corresponding  times  for  central 
differences . 

This  work  will  appear  in  a  Ph.D.  dissertation  by  E.  W.  Fieri 
(in  preparation),  and  the  two-dimensional  version  of  the  finite 
fluid  element  method  will  also  appear  there. 
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7.  Vortex  motion 

In  some  problems  of  fluid  mechanics,  there  is  a  boundary 
Integral  method  that  can  be  used  to  formulate  efficient  numerical 
methods.  A  particular  example  is  the  flow  induced  by  a  two- 
dimensional  region  of  constant  vorticity,  i.e. 

^  ^  =  a)  in  R  bounded  by  C 

=  0  outside  R 


The  solution  of  this  problem  for  which  u  — *•  0  as  |x|  — ►  <*>  is 


u(X,Y,t  )  =  f-  /  —X  ~  Yg  dxdy 

R  lx  -  xr 


v(x’Y't)  =  -^lTrf?dxdy 


and  this,  in  turn,  can  be  transformed  to  the  line  integral 


u(X,t)  =  "  2?^  £nl*  ~  X|d*  • 

C 

For  two-dimensional,  incompressible  flow  du/dt  =  0  on  particle 
paths,  so  the  motion  of  the  boundary  of  R  is  governed  by 


X  =  -  ^  f  in|x  -  X|d(x-X)  , 
C 


X  on  C. 


The  shift  of  the  integration  variable  from  x  to  x-X  removes  the 
logarithmic  singularity,  for  after  an  integration  by  parts  on 
the  closed  curve. 
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The  normal  component  of  the  velocity  of  the  boundary  point  X  is 


A  • 

N*X 


=  —  f 

2  it  J 


N • ( x  -  X ) ( x  -  X)  •  dx 


lx  -  X 


A  • 

and  the  contribution  to  N*X  at  x=X  is  zero. 

There  are  a  number  of  ways  to  perform  the  numerical  integration 
the  one  we  have  set  up  fits  a  periodic,  cubic  spline  to  the 
curve  x(  e),  i.e. 

N  u 

*(  9)  -  y  a1B1(  0)  +  0(h*) 


where  Bj^  is  the  cubic  B-spline  and  h  is  [x^  -  x1+1  I-  The 
derivative 

N  o 

d_x  =  y  ^B^ejde  +0(h-5de) 


gives  the  tangent  vector,  the  normal  is  obtained  by  rotation, 
and  the  normal  component  of  X  is  calculated  in  a  consistent 
approximation  by  Simpson's  rule. 

With  one  iteration  of  central  differences  for  the  time 
integration  (Huen's  method)  the  algorithm  performs  well,  but 
we  are  not  yet  satisfied  with  it.  Further  work  on  this  subject 
is  planned. 


8.  Graphics 

Visualization  of  two  and  three  dimensional  flows  creates 
special  problems  and  a  need  for  special  software  that  isn't 
generally  found  in  standard  packages.  Graphics  functions  we 
developed  for  our  own  purposes  include: 

1)  Standard  plotting  of  rough  graphs  at  a  terminal  or 
better  quality  graphs  onsn  X-Y  plotter 

2)  Rough  representation  of  direction  fields  and  contour 
maps  at  a  terminal. 

3)  Contour  maps  on  X-Y  plotters. 

4)  The  linear  algebra  (affine  transformations  and 
catenations)  of  two-dimensional  directed  arcs. 

5)  The  linear  algebra  of  space  curves,  including  orthographic 
and  stereographic  projections  and  binocular  pairs  of 
projections. 

Item  4)  was  of  such  a  general  utility  that  it  has  now  been 
included  in  the  public  APL  software  library  at  Brown  University. 
Items  3)  and  5)  will  probably  be  included  in  the  public  domain, 
as  well,  but  in  any  case,  listings  of  any  of  our  graphics 
software  functions  are  available  on  request. 

The  most  challenging  graphics  problem  was  the  generation  of 
contour  maps.  By  comparison  with  other  contour  map  algorithms 
we  have  seen,  our  approach  is  somewhat  different,  and  it  appears 
to  be  considerably  more  efficient.  A  description  of  the  contour 
map  algorithm  is  included  here  in  appendix  P. 


Appendix  A 

Abstract 

This  paper  contains  the  formulation  of  a  numerical  method 
to  treat  one,  two  or  three  dimensional  unsteady  flow.  The 
method  is  closely  related  to  the  PIC  method  and  finite  element 
methods.  The  elements  (finite  fluid  elements)  are  localized 
distributions  of  mass,  characterized  by  mass,  radius  and  position 
in  space.  They  move  according  to  interactions  between  neighboring 
elements  that  are  derived  from  the  fluid  equations.  Rough 
estimates  of  the  dependence  of  the  accuracy  of  the  approximation 
on  particle  radii  and  the  number  of  neighbors  retained  are 


calculated. 
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1 .  Introduction 

This  paper  is  devoted  to  the  formulation  of  a  numerical  method 
to  treat  one,  two,  or  three-dimensional  unsteady  flow  of  a 
barotropic,  Stokes  fluid.  The  generalisation  of  the  method  to 
treat  an  arbitrary  compressible  fluid  will  be  deferred  until  the 
present  model  has  been  implemented,  and  its  worth  as  a  practical 
algorithm,  appraised. 

The  method  to  be  developed  here  is  closely  related  to  the 
PIC  method  (Particle  In  Cell,  Refs.  l-3)»  but  in  fact,  the 

two  essential  features  of  that  method  will  both  be  treated  somewhat 
differently.  The  particles  will  not  be  treated  as  point-masses; 
instead  they  will  be  taken  to  be  localized  distributions  of  mass 
of  finite  extent,  and  as  we  shall  see,  they  will  be  allowed  to 
overlap.  The  addition  of  another  parameter,  an  effective  radius, 
for  each  particle  allows  another  degree  of  freedom  for  optimiza¬ 
tion,  and  thus  an  inherently  more  accurate  representation  of  a 
continuous  density  field  by  a  finite  set  of  localized  particles. 

The  equations  of  motion  of  the  particles  are  derived  directly 
from  the  equations  of  fluid  mechanics,  and  the  constitutive 
relations  of  the  fluid  are  modelled  by  providing  appropriate 
interactions  between  particles.  Here,  as  in  the  PIC  method, 
closure  of  the  dynamical  system  of  particle  motions  is  effected  by 
fitting  mass  and  momentum  density  fields  to  provide  a  continuum 
approximation  of  the  motion  of  the  particles.  The  system  is 
then  closed  by  adopting  the  constitutive  relations  of  the  fluid 
that  is  to  be  described.  Instead  of  the  fixed  cells  of  the  PIC 


method,  cells  that  move  with  the  particles  will  be  employed. 


Mass  and  momentum  density  in  the  neighborhood  of  each  particle 
will  be  defined  by  weighted  averages,  over  the  particle  and  its 
neighbors,  of  the  corresponding  particle  attributes.  This  means 
that  the  implementation  of  the  method  for  two-  and  three-dimensional 
flows  will  require  the  tabulation  of  a  continually  updated  list  of 
the  near  neighbors  of  each  particle.  It  does  not  appear  that  this 
is  a  burden,  however,  since  the  major  contribution  to  the  weighted 
averages  comes  from  the  six  or  twelve  nearest  neighbors  in  two  or 
three  dimensions,  respectively,  and  the  updating  can  easily  be 
accomplished  by  periodically  checking  the  positions  of  second 
neighbors,  which  are  easily  assembled  from  the  neighbors’  neighbors. 

The  method  is  also  closely  related  to  the  finite  element 
methods  that  have  been  used  in  continuum  mechanics  (Refs.  ^,5). 
However,  two  departures  from  the  usual  scheme  of  things  will  be 
introduced:  In  the  first  place;  the  contiguous,  nonoverlapping 

elements,  typically  triangles  or  quadrilaterals  in  two  dimensions 
and  tetrahedra  and  various  prisms  in  three,  will  be  abandoned  in 
favor  of  simpler  elements,  the  finite  particles,  that  fill  space 
in  an  additive  way  rather  than  by  mutual  exclusion;  i.e.  they 
overlap.  Thus  the  element  parameters,  nodal  values,  coefficients 
of  shape  functions,  etc.,  are  replaced  by  particle  attributes  and, 
as  we  shall  see,  cell  attributes.  The  basic  similarity  between  this 
and  other  finite  element  methods  is  in  the  determination  of  the 
particle  and  cell  attributes  by  Galerkin's  method  where  the 
integrated  square  of  the  error  of  the  finite  approximation  is 


minimized . 


The  second  departure  from  other  finite 
element  methods  is  that  this  one  is  entirely  a  Lagrangian 
formulation:  the  velocity  of  each  particle  is  among  the  attributes 

to  be  defined  by  minimization  of  the  error  of  the  approximation. 

The  velocities  are  then  to  be  integrated  to  find  the  positions  of 
the  particles,  thus  giving  a  formulation  that  closely  mimics  the 
fluid  elements  on  which  fluid  mechanics  is  founded. 

With  the  overlapping  of  particles  and  the  use  of  Galerkin's 
method  to  govern  the  coefficients  of  an  additive  covering  of  space, 
the  present  method  bears  at  least  a  superficial  resemblance  to 
spectral  methods  that  have  been  used  in  numerical  fluid  mechanics 
(Ref.  6).  In  the  spectral  methods  the  'elements' 

are  trigonometric  functions  and/or  members  of  various  sets  of 
orthogonal  polynomials,  none  of  which  are  localized  in  space. 

The  entity  that  corresponds  to  an  inner  product  (here  denoted  by 
<f^f^>)  is  the  integral  over  all  space  of  the  product  of  the  ith 
and  jth  mass  distributions.  The  normalization  is  different, 
however;  it  is  <f^>  -  1  here.  The  matrix  <f^fj>  changes  with  time, 
and  even  the  location  of  the  largest  off-diagonal  elements  changes 
in  a  shear  flow;  thus  its  inversion  has  to  be  done  numerically. 
Because  of  the  localization  of  the  particles,  however,  the  major 
off-diagonal  contributions  to  the  i^  row  come  from  relatively  few 
near  neighbors  of  the  i^*"1  particle. 


2 .  Fluid  Equations 


Ait 


The  fluid  model,  a  barotropic,  Stokes  fluid  includes  some 
features  of  compressible  flow.  Specifically,  the  gradient  of 
the  pressure  is  retained  in  the  momentum  equation,  but  the 
thermodynamics  is  simplified  by  taking  the  pressure  to  be  a 
function  of  density  alone.  The  Stokes  approximation  defines  stress 
in  terms  of  a  single  viscosity  coefficient,  also  a  function  of 
density,  and  the  theory  is  closed  without  the  need  for  an  energy 
equation.  Enough  of  the  complexity  of  a  compressible  fluid  is 
retained,  so  that  the  complementary  functions  of  the  particles 
and  the  cells  in  the  finite  fluid  element  model  can  be  fully 
appreciated. 

The  equations  that  have  a  tensorial  character  will  be  given 
twice,  first  in  the  cartesian  tensor  notation  and  then  in  an 
abbreviated  notation  that  will  be  used  in  later  sections  where 
the  subscripts  will  denote  particle  and/or  cell  number.  The 
cartesian  tensor  subscripts  take  on  one,  two  or  three  values  in 
as  many  space  dimensions,  summation  over  repeated  indices  is 
implied,  and  6^  is  one  if  i=j ,  zero  otherwise. 

Conservation  of  mass: 


Ip.  +  9 _ 

at  axj 


(PUj) 


o 


(i) 


Dj.  +  V '  ( pu )  =  0 
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Momentum : 


(2) 


ft  ((JU1>  +  Hr  <CUJU1)  ■  cgi  +  hr  °J1 


(pu)t  +  V(puu)  =  pg  +  Vo 

where  g  is  a  prescribed  external  force  per  unit  mass,  and  the 
stress  (with  Stokes'  approximation)  is 

3u  3u  3u, 

°lj  =  u(JZJ  +  3Xi  “  3(3^)6ij}  ‘  p6ij 


(3) 


a  =  2yVu  -(p  +  yV'u)I. 


Finally,  for  the  barotropic  model: 


(M 


P  =  P(p) 

u  =  u(p) 


where  the  functions,  p  and  u,  are  prescribed. 


3 •  Particle  attributes 

Just  as  in  other  finite  element  methods,  the  individual 
elements,  particles  in  this  case,  can  be  assigned  an  arbitrarily 
complex  internal  structure.  This  leads  to  the  usual  trade-off: 
the  accuracy  of  the  approximation  can  be  improved  either  by 
increasing  the  number  of  elements  per  unit  volume  or  by  augmenting 
the  internal  complexity  of  the  elements.  The  choice  is  by  no 
means  a  trivial  one,  for  though  the  more  complex  elements 
necessitate  more  equations  to  govern  the  values  of  their  several 
attributes,  the  greater  separation  of  elements  tends  to  allow 
larger  time-increments  in  the  numerical  integration  of  the 
dynamical  system.  The  question  of  improvement  of  the  approxima¬ 
tion  will  be  deferred  here,  and  the  particles  will  be  assigned  a 
relatively  simple  internal  structure. 

The  particle  approximation  consists  in  the  replacement  of  the 
density  field,  p(x,t),  by  a  set  of  N  localized  distributions  of 
mass ,  i . e . 

(1)  p  ~  Em1(t)f1(x,t) 

where  the  sum  is  from  i=l  to  N.  The  identification  of  m^  as  the 
total  mass  of  the  i^  particle  is  effected  by  the  conditions  on 
the  distribution  functions, 

(2) .  <f±>  =  1  i=l, .  .  .N, 

where  the  angular  brackets  denote  the  integral  over  all  space. 
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Since  the  distribution  functions  represent  the  location  of  mass 
within  the  particles,  the  location  of  the  ib  particle  is  the 
mean  value. 


(3)  X±(t)  =  <xfi>, 

4*  V> 

and  the  effective  radius  of  the  iL  particle  is  defined  in  terms 
of  the  variance  as 

(4)  r^t)  =  <|x-X1|2f1>1/2  . 


If  the  attributes,  mass,  position  and  effective  radius,  are 
accepted  as  a  sufficient  description  of  the  particles,  then  a 
particularly  convenient  choice  of  distribution  functions  is  the 
normal  distribution. 


(5)  =  ^u/2 

where  the  value  of  n  is  1,  2 


- | x-X± | 2/a| 
e 

or  3  in  n  space  dimensions. 


and 


(6)  cr^t)  =  J2-  r1(t). 

It  may  be  noted  that  the  first  few  moments,  eqs.(2,3»*0,  are  far 
from  sufficient  information  to  determine  the  distribution  functions 
the  arbitrary  and  convenient  choice  of  normal  distributions  relates 
higher  moments  to  the  radius  in  ways  that  are  particularly  easy  to 
evaluate.  A  greater  convenience,  however,  is  the  ease  with  which 
several  kinds  of  matrix  elements  that  will  be  introduced  in  due 
course  can  be  evaluated  explicitly.  Thus,  for  example,* 

if 


See  Sec.  6,  Matrix  elements. 


With  the  normal  distributions,  eq.(5),  the  mass  of  the  ith 
particle  is  distributed  symmetrically  about  Xi(t),  and  the 
approximation  to  the  fluid  momentum  field  that  is  consistent 
with  eq. (1)  is. 


(8)  pu  ~  Em^X^f^ 

where  u(x,t)  is  the  fluid  velocity  field.  Clearly,  X^  and  X^ 
cannot  be  specified  independently,  so  to  preserve  the  character 
of  the  initial  value  problem  for  the  fluid,  the  parameters, 
m^  and  X^  are  chosen  to  minimize  the  error  of  the  approxima¬ 
tions  of  eqs.(l  &  8).  For  the  particle  approximation  the  error 
will  be  taken  to  be  the  unbiased  r.m.s.  error, 


(9) 


E  =  ~  <(p-Zmifi)2  +  |pu-Em1Xifi|2>1/2  , 

U* 


where  U#  is  a  characteristic  fluid  velocity  that  has  been  intro¬ 
duced  to  render  eq.(9)  dimensionally  consistent,  p#  is  a  character 
istic  fluid  density,  and  [E  ]  is  volume.  The  choice  of  eq.(9) 

is  somewhat  of  a  compromise,  for  though  an  unbiased  relative  error 
e.g. 


(10) 


(p-Em^)2  |pu-Zm1X1fi  |2 

<  —  -f  —  >  , 

P*'  I  pu  |  2 


might  provide  a  better  representation  of  the  fields,  the  theory 
that  follows  therefrom  is  much  more  difficult. 


A9 


The  conditions  to  minimize  E  are: 


(ID 


p2e  . 


3E 


p^EpL  =  <Pf  >  -  £m . <f . f . >  +  X  •  =  0 

3m1  i  J  i  J  mi  i  3X^ 

2  2 


<i2)  ^Ef?-  -  <pufy  -  £Vj  “W  * 0 


'i  3X 


(13)  ^  +  ^  xl-<f’c1<pu-!:"’jXJfj)> 


1  - 

2 

* 


=  0 


where 


(14) 


3f 


ai  '  9ai  o\ 


.  O  O  nC?4 

1  =  ^  (Ix-xj2  -  ~2i)fi  • 


Note  that  eq.(12)  implies  the  vanishing  of  the  corresponding  term 
of  eq. (11) ,  and  eqs.(ll  and  12)  allow  the  replacement  of  f^  by 
|x-Xi|2f1  in  eq . (13 ) • 

Now  it  may  be  noted  that  if  the  conditions  of  equation  (13) 
are  retained,  the  result  is  a  theory  in  which  the  cr's  are  large , 
the  particles'  are  not  localized,  and  near  and  far  neighbors  are 
of  equal  importance  in  equations  (11)  and  (12).  Accordingly, 
equation  (13)  will  be  dropped  altogether,  and  in  §7  estimates  will 
be  given  for  choices  of  a’s  that  minimize  the  error  of  truncated 
approximations . 
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4 .  Particle  dynamics 

Consider  first  eqs . ( 3 . 11 , 3 • 12 ) : 

(1)  Enij  <f  j  >  =  <pf1> 

(2)  EmjXj<f^fj>  =  <pufi>  . 

In  the  description  of  unsteady  flow,  the  parameters,  m^, 
and  o^,  are  to  vary  with  time,  but  in  such  manner  that  eqs. (1,2) 
are  preserved.  Equations  (1,2)  can  thus  be  regarded  as  constraints 
on  the  particle  dynamics;  to  see  their  effects,  it  suffices  to 
differentiate  them  to  obtain 

(3)  (EmJ<f1fj>)  =  <pf1>*  =  <pfi+fiPt> 

(4)  (EmJXj<fifJ>)'  =  <puf1>'  =  <pufi+fi(pu)fc>  . 

The  subscript,  i,  runs  from  1  to  N  in  eqs. (3,4),  which  can  be 
thought  of  as  the  mass  and  momentum  transport  laws  of  the  particle 
dynamics.  The  right-hand  sides  of  eqs. (3,4)  contain  both  particle 
attributes  and  the  fluid  density  and  momentum  fields;  they  are  the 
particle  interactions  that  are  to  be  determined  presently  by  the 
introduction  of  cells. 

From  the  fluid  eqs . (2 . 1 , 2 . 2 )  and  eq.(3-l4)  it  follows  that 

•  2 

2a,  p  no, 

(5)  <pf1>  =  — |  <p(|x-Xir  -  "2i)fi> 

ai 

-  <pVf1>  •X1  -  <fiV*pu> 


All 


2 

2a,  ?  no. 

(6)  <puf1>'  =  — j  <pu(|x-Xi|  - 

°i 

-  <puVf^>,Xi  -  <f^7*puu>  +  <f^(pg+V-o)> 

Clearly,  eqs. (3, 4, 5, 6)  are  incomplete,  even  if  the  mass  and 
momentum  densities  were  known;  there  are  only  2N  relations  for 
the  time  derivatives  of  3N  particle  attributes.  The  remaining 
N  equations  can  be  obtained  by  differentiating  the  estimates  of 
o^t)  that  will  be  given  in  §7,  o£  eqs.  (3, 4, 5, 6)  can  be  iterated 
with  =  0  in  the  leading  approximation. 


Some  estimates  of  the  relative  importance  of  various  terms 
in  eqs. (5,6)  can  be  made  as  follows:  Let  it  be  supposed  that  the 
Taylor  series  of  the  mass  and  momentum  densities  were  known,  i.e. 

(7)  P  =  Pi+(x-X1)  •  (Vp)1  +  |(x-X1 )  (x-Xj^ )  :  (VVp^  +  ... 

(8)  pu  =  (pu)1  +  (x-X1)  •  (Vpu)i  +  |-(x-Xi)  (x-Xj^)  :  (VVpu)i  +... 

where  the  subscripted  quantities  are  evaluated  at  x  =  X^.  The 

substitution  of  the  ittl  Taylor  series  in  the  ith  of  eqs.  (5,6), 

an  integration  by  parts  of  the  matrix  elements  that  contain  Vfi, 

and  the  results  for  normal  distributions,* 

2 

P  no, 

(9)  <|x-Xi|2f1>  = 

- - 


See  Sec.  6,  Matrix  elements. 
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2 

(10)  <  I  x-Xi  I  V±>  =  (J  +  g-)oJ  , 

give  the  results, 

(11)  <pf1>‘  =  (~  a1a1(Ap)i+(7p)i-Xi-(\7-pu)1)(l+O(0^)) 

(12)  <pufi>‘  =  (|  a1ai(Apu)i+(Vpu)1*X1-(V* puu)1 

+Pigi+(V-a)i)(l+0(o^)). 

The  coefficients  of  the  relative  rates  of  change  are  of 

the  order,  0(a.),  and  thus  as  the  number  of  particles  per  unit 
volume  is  increased,  the  error  introduced  by  neglecting  those 
terms  of  eqs. (11,12)  vanishes  as 

(13)  o\  =  0((V/N)2/n)  -  0  . 

f 

In  the  same  limit,  approaches  as  N/V  becomes  large,  and 
the  major  contributions  to  the  particle  dynamics  are: 

(1*0  (ImJ<f1fJ>)  +  Pi(V*u)1  =  0 

(15)  (LmjXj<fif  j>)  +  (pu)i(V-u)i  =  p^  +  (V-o)1 

where,  once  again,  the  subscripts  on  the  fields  and  their  gradients 
indicate  evaluation  at  x  =  X^(t). 


kJ _ _ 


From  the  results  of  the  previous  section  it  can  be  seen  that 
the  least  information  about  the  fluid  mass  and  momentum  densities 
that  is  needed  to  complete  the  particle  dynamics  approximately 
includes  the  values  at  x  =  X^,  and  (pu)^,  the  first  derivatives, 
(Vp)i  and  (Vpu)^,  for  the  evaluation  of  (V*u)^,  and  the  second 
derivatives,  (VVp)^  and  (VVpu)..,  for  the  evaluation  of  (V-a).. 

For  that  purpose  then,  cells  are  introduced  to  provide  local 
polynomial  approximations  where 

(1)  P  ~  Pl  +  (x-X1)-(Vp)i  +  |(x-Xi)(x-Xi):(VVp)i 

(2)  pu  ~  (pu)1+(x-X1)  •  (Vpu)1  +  i(x-Xi)  (x-X.^)  :  (VVpu^  . 

The  subscripted  quantities  in  eqs.(l,2)  are  not  the  fields  and 
their  derivatives  evaluated  at  x  =  X^,  as  in  the  previous  section; 
rather,  they  are  the  coefficients  of  quadratic  approximations  of 
mass  and  momentum  density  fields  that  would  be  obtained  if  the 
system  of  particles  were  treated,  after  the  manner  of  the  kinetic 
theory  of  gases,  in  a  continuum  approximation.  The  process  is 
the  reverse  of  kinetic  theory,  however,  since  it  is  the  particle 
interactions  that  are  unknown  here,  and  the  mass  and  momentum 
densities  are  to  be  found  in  order  to  adjust  the  particle  inter¬ 
actions  to  suit  the  thermodynamic  properties  of  the  fluid  at  hand. 

The  local  quadratic  approximations  for  the  mass  and  momentum 
densities  are  determined  by  minimization  of  the  biased  r.m.s.  error. 


(3) 


Al4 


_  1 


P* 


;Fi[(p1+(x-Xi) * (Vp)1+  |(x-X1)(x-X1):(7Vp)1-ZmJfJ )2+ 


+  ( (pu)1+(x-Xi)  •  (Vpu)i+  ^-(x-X^  )  (x-X^ ) :  ( VVpu 

u# 

-  E”jijfj)2]>1/2 

where  F^(x,t)  is  a  localized  weight  function  that  is  normalized, 

(4)  <Ft>  =  1, 
located  at  X  (t), 

(5)  Xi(t)  =  <xFi>  , 


and  has  an  effective  radius, 

(6)  R1(t)  =  <|x-X1|2F1>1/2. 


Again,  for  the  ease  with  which  matrix  elements  can  be 
evaluated  explicitly,  the  weight  functions  are  taken  to  be  the 
normal  distributions. 


(7) 


F, (x, t )  = 


“Ix-Xj^  |  2/E2 


(  ttz|  ) n/2 


where  n  is  the  number  of  space  dimensions  and 


(8) 


Z1  ( t ) 


1 


-  R, 
n  i 


(t). 


This  time,  the  quantities  that  enter  in  the  conditions  to 
minimize  are  the  subscripted  field  quantities  and  E  ^ ;  the 
conditions  are: 


(9) 


+  i  <(x-X1)(x-X1)F1>:(V7p)1  =  EmJ<FifJ> 


(10) 


( pu )  ^  ^  <(x-Xi)  (x-X1  )Fi>  :  (VVpu)1  =  Zrrij  Xj  <F±  f  ^  > 

(11)  <(x-X1)(x-X1)F1>-(Vp)1  =  EmJ<(x-X1)F1fj> 

(12)  <(x-Xi)(x-Xi)Fi>-(Vpu)1  =  EmjXj<(x-Xi)FifJ> 

(13)  <(x-Xi)(x-Xi)Fi>Pi  +  |<(x-X1)(x-Xi)(x-Xi)(x-Xi)Fi>:  (Wp)i 

=  EmJ<(x-Xi) (x-Xi)F1fj> 

(14)  <(x-X1)(x-Xi)Fi>(pu)i  +  |<(x-Xi)(x-Xi)(x-X.)(x-X1)Fi>:(VVpu) 

"  ^jXj<(x-Xi)(x-X1)Fifj> 

and 

(15)  <FJ.^[(p1+(x-Xi)  •  (Vp)1+  |(x-X1)(x-Xi):(Wp)1-Emjfj.)2  + 

+  ^  ( (pu)1  +  (x-X1)  •  (Vpu)1+  |-(  x-X^ )  (x-X^ )  :  (VVpu)^ 

-  Em.X.f, )2]>  =  0 
j  J  1 

where 

9Fl  p  ?  nEi 

0.6)  Fj.^  =  -gy-  =  p-  (|x-X1|  2~)pi' 

In  eqs.(9  to  14)  the  terms  containing  odd  moments  of  F^  have 
already  been  dropped;  in  the  next  section,  further  simplifications 
that  follow  from  the  choice  of  normal  distributions  will  be 
derived . 

The  final  condi¬ 
tion,  eq.(15)  which  governs  the  choice  of  cell  radii  to  minimize 
the  error  of  the  approximation,  implies  the  inappropriate  result, 
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=  0-  Accordingly,  it  will  be  dropped,  along  with  the  corre¬ 
sponding  condition  of  eq.(3-13),  and  other  estimates  of  both 

and  will  be  given  in  the  section  on  particle  and  cell  radii. 

Given  the  quadratic  approximations  for  the  mass  and  momentum 

densities,  the  corresponding  fluid  velocity  is 

(17)  u  =  ui  +  (x-Xi)  •  (Vu)^  +  i-(x-X1)(x-X1):(VVu)i  +  0  (  |  x-X±  |  3  ) 
where 

(18)  u.^  =  (pu;i/pi 

(19)  (Vu)±  =  ((Vpu)1  -  u1(Vp)1)/pi 
and 

(20)  (VVu)1  =  ((VVpu)1-ui(77p)1-(Vu)1(Vp)i-(Vp)i(Vu)1)/pi. 

The  thermodynamic  quantities  that  appear  in  Vo  are 

(21)  p  ~  p  (  p1  )+p  ’  (  p^  (x-X1)  •  (Vp)i  + 

+  |(x-X1)(x-X1):(p'(Pi)(VVp)1+p"(P1)(Vp)1(Vp  )i) 

and 

(22)  y  ~  y(pi)  +  UT ( p± ) • (Vp)±  + 

+  i(x-X1)(x-Xi):  (u1  (p1)(VVp)1+v"(P1)(Vp)1(Vp)i). 
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6.  Matrix  elements 


The  choice  of  the  normal  distributions, 
1 


(1)  fl  =  („,2)n/2  e 


-  I  x-X± | 2/o\ 


(2)  Fi 


(tte|  )n/2 


leads  to  numerous  simplifications  of  the  several  matrix  elements 
that  have  been  introduced.  Generally  speaking,  the  simplification 
of  the  tensorial  character  of  the  matrix  elements  is  a  consequence 
of  the  choice  of  functions  of  |x-X^|  ,  and  the  relative  ease  with 
which  they  can  be  evaluated  is  a  consequence  of  the  specific  choice. 

Consider  first  the  right-hand  side  of  eqs.(5.9  to  S-l2*)-’ 

The  integral  <F^fj>  is  of  the  same  form  as 


(3) 


<fiV 


-|x-Xi 


2/  2 

/oi 


- | x-X4 


,  2  2  2 .n/2 
(n 


2/  2 

/0J 


dx 


and  th^  integral  <F^ Tj  f'1<>  that  will  appear  in  the  next  section. 
The  integral  of  the  product  of  any  number  of  normal  distributions 
can  be  evaluated  by  introducing  the  ’center  of  mass’, 

X,  . 

(A)  X  =  (I  -|)7e  ±_  • 

°i  ai 

The  exponent  in  the  integrand  is  then 


|x-X|2 


|X-X±|2 
- 2 - 


+2(x-X) 


(X-X  ) 

Z  - . 


o 


i 


The  last  term  vanishes  by  the  choice  of  X,  the  second  term  gives 
the  dependence  on  the  coordinates  X^,  and  the  integral  of  the  first 


l8 


term  gives  the  normalization.  The  result  for  the  product  of 
two  can  be  rearranged  to  give  eq.(3-7)  and 

-  |  X.-X  | 2/( E^+a?) 

(6)  <Fif1>  =  - 2-  g  n/2  e  J  • 

1  J  (it  (Ej+oJ))n/<f 

Likewise  the  result  for  the  product  of  three  normal  distributions 
can  be  rearranged  in  the  form 


<FifjV 


,  2, _2,  T~T7~2~277n72 
(w  (I1(oJ+<Ik)tojok)) 


Ei|xJ-xk|S2|vxi|,;+ok|xj-xi 


2,2 


Ei(0j+°k)+°jak 


The  rearrangement  of  the  general  result  is  somewhat  more  tedious 
than  the  direct  evaluation  of  eq.(7)  by  completion  of  the  square 
in  the  exponent. 

Next,  in  order  of  appearance,  is 


(8)  <(x-X1)F1f  >  =  <F1f1 


'  J-2  (xrxi)  <FiV  > 


and  then 


Zi+0j 


(9)  <(x-Xi)(x-Xi)FifJ>  =  ~  (<(x-XJ|  )F<f<>Y  +  I<FJ,>) 


i  i  j  X, 


=  (  -3-  ?  ?  (x,-x,  ) (X, -X,  )  + 

(Ei+aj)  J  J 

y2  2 
E .  a . 

+  -  V  g-  I)  <F .  f .  > 

2 (E?+o?)  1  J 
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where  I  is  the  unit  dyad  (equal  <5^  in  cartesian  tensor  notation). 

On  the  left-hand  sides  of  eqs.(5.9  to  5.1*0  the  matrix 
elements  are: 

(10)  <(x-Xi)(x-X1)Fi>  =  -§  1 
and 

(11)  <(x-X1)(x-X1)(x-Xi)(x-X1)Fi>  =  J  ijs 

where  S  is  the  symmetric  sum  of  outer  products  of  unit  dyads 
that  has  the  cartesian  tensor  representation. 


’i  j 


6ij6kJ. 


6i*6jk 


The  results  (10,11,12)  follow  by  the  substitutions  of  (x-X^  and 

(x-X.^)  (x-X^)  (x-X^)  for  0 ( x )  in  the  identity, 

2  2 
-l  1 

(13)  <(x-X1)F1G>  =  <VFiG>  =  <Fj,VG>  . 


Contractions  of  (10,11,12)  then  give  the  results  (4. 9, *1.10). 

This  completes  the  evaluation  of  matrix  elements,  except  for 
moments  of  fj f^>  that  can  be  evaluated  after  the  manner  employed 
for  eqs.(8,9).  With  the  simplification  that  follows  from  the 
choice  of  normal  distributions,  eqs  .  ( 5  •  9-5 • 14 )  now  are: 

(14)  p±  +  J.  L^(Ap)i  =  Imj  <F1fJ> 

X  —X 

(15)  (7p),  -  SEmj  4-4  <F1fJ> 
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(16)  PjI  +  \  z|( ( Ap)iI+2(77p)1)  = 


Z?(X  -X,)(X.-X.)  a?I 
=  2Em.  (  -  — *2 —  +  -  % 


(Zi+aj}‘ 


-)  <F,  f  .  >  . 


2(E^+oJ) 


From  eq.(l^J)  and  the  trace  of  eq.(l6)  it  follows  that 


m,  X.-X, 


(17)  (Ap)i  -  Hi  2J2  (  ^  2  J)  <Fifj> 


Ei+0j  Ei+0j 


Ef  | X , -X . 

(18)  p  =  Em.  (1  - 

J  Ei+aj  Ei+°J 


The  corresponding  relations  between  coefficients  for  pu  are 

identical  except  that  m.  is  everywhere  replaced  by  m.X.. 

J  J  J 


2j 


7 •  Particle  and  cell  radii 

At  the  end  of  §4  it  was  mentioned  that  the  attempt  to  minimize 
the  unbiased  error  by  choice  of  o's  leads  to  a  nonlocalized  theory 
in  which  the  ’particle'  radii  are  large.  The  result  can  be  seen 
by  considering  the  1-dimensional ,  uniform,  periodic  case  where 
P=1 ,  -oo  <  x  <  »,  the  particles  are  equally  spaced  at  =  i  with 
equal  masses  m^  =  m,  and  =  a.  Then 

(1)  E2  =  <(l-mEf  )2> 

and  the  relation  between  m  and  a  is 

oo  2  2 

(2)  <f  >  =  1  =  mE<f . f , >  =  -^=-  Z  e-k  /2a  . 

1  “  /2tt o 

Prom  equations  (1,2)  it  follows  (as  in  the  derivation  of  Bessel's 
inequality  for  orthogonal  functions)  that 

(3)  E2  =  <1>  -  2m  E<f±>  +  m2  ZZ  <fAfj> 

=  <1>  -  m2  EE<fifJ> 

=  <l-m>  . 

Since  the  particles  are  equally  spaced  the  error  per  particle  is 

o  00  |2 

(*1)  e  =  1-m  =  1  -  /2tt  a/  Z  e  ' d0  . 

—  OO 

o 

Figure  1  is  a  graph  of  e  for  0  <  o  <  1,  showing  its  decline  to 

2 

zero  as  a  ■*  ®  :  the  notable  feature  is  the  rapid  approach  of  e  to 
negligibly  small  values.  A  similar  phenomenon  takes  place  in  two 


and  three  dimensions. 


Given  that  a's  will  be  relatively  small,  the  next  question 
is  the  truncation  error  associated  with  the  approximate  evaluation 
of  m’s  from 

(5)  <pf,>  ~  £  m.<f . f . > 

where  the  sum  is  over  j=i  and  the  N  nearest  neighbors  of  the  1th 
particle.  In  the  uniform  periodic  case  the  m's  and  o's  are  equal, 
p  is  1,  and  the  particles  are  at  the  integers  in  one  dimension, 
in  a  hexagonal  lattice  in  two,  and  in  one  of  the  close-packed 
lattices  in  three.  In  any  case 

(6)  1  =  mM(o)  £  <f  f,>  , 

<  N  °  J 

(7)  E2  =  <1>  -  2mN  Z<f1>  +  m2  I  E  <fif1> 

OO  00  00  1 

=  <1>  -  £  mN  +E  £  m2  <f  f,>, 

OO  OO  >  f\J  ^ 

and  the  error  per  particle  is 

(8)  =  Tn-„,H(l-mN  ^  <f0fj>) 

where  xn  is  the  volume  per  particle  in  the  n-dimensional ,  close 
packed  lattice,  i.e. 

(9)  x  =  d ,  ^d2,  -  d3  for  n  =  l,2,3 

"  2  /2 

where  d  is  the  separation  of  nearest  neighbors.  Figures  2,3  and 

p 

4  are  graphs  of  e  for  truncations  that  include  first,  second  and 
third  neighbors,  for  various  ranges  of  a.  The  separation  of  nearest 
neighbors  is  set  equal  to  one  in  all  cases,  and  the  three-dimensional, 
close-packed  lattice  is  face-centered  cubic. 
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A  different  measure  of  the  accuracy  of  the  truncated 

approximations  is  given  in  Figures  5,  6  and  7,  which  are  graphs 

of  Tn-MN.  The  value  of  a  at  which  =  xn  is  the  one  for  which 

mass  is  conserved,  i.e.  <p>  =  Em^ .  It  should  be  noted,  however, 

that  the  value  of  a  at  which  mass  is  conserved  is  systematically 

2 

less  than  that  where  e  is  minimized.  By  direct  calculation,  the 
values  of  a  for  mass  conservation  with  the  nearest  two,  six  and 
twelve  neighbors  included  are: 
a 

(10)  s  0.572,  0.498,  0.444  for  n=l,2,3- 

The  corresponding  errors  per  unit  volume  are: 

2 

e 

(11)  ^  dn  s  0.000311,  0.00887,  0.0258. 

n 

In  the  case  where  p  is  not  constant  and  the  particles  are 
not  in  a  periodic  array,  a  simple  estimate  of  cr^t)  is  equation  (10) 
with 

(12)  d  (t)  -  Z  | X .-X . |/N 

1  <N  J 

where  the  sum  is  over  nearest  neighbors.  Many  other  ways  to 
estimate  can  be  found,  however,  and  the  adoption  of  (12)  may 
be  considered  as  provisional.  In  any  case,  as  the  number  of 
particles  per  unit  volume  becomes  large,  -*•  0  along  with 

d^,  and  If  p  is  differentiable,  then  equation  (11)  provides  an 
estimate  of  a  local  relative  error 

(P-Em.f , )2 

(13)  < - J-*3 —  >1/<1>i 

p 
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where  <1>^  is  the  volume  of  a  neighborhood  of  X^  that  has  a 
diameter  that  is  large,  compared  to  d^,  and  is  small,  compared 
to  the  scale  length  p/|Vp|. 

There  remains  now  the  question  of  the  error  in  the  cells. 
The  uniform,  periodic  case  gives  rise  to  double  sums  that  are 
relatively  easy  to  evaluate  in  one  dimension,  rather  more 
difficult  in  two  and  three.  The  cell  error  is 


<Fo(pirmIfj  }  > 


where 


-|X  |2/(Z2+o2) 

m  £<F  f . >  =  - — —pr  l  e  J 

<N  0  J  ( 7T  (X2+a2))n/2  <N 


and  again  the  sum  is  over  the  N  nearest  neighbors.  From  (14,  15) 
it  follows  that 

(16)  E2q  =  m2E  £<F  f.f .>  -  PN(pNv2m  E  <F  f  >), 
and  in  the  one  dimensional  case 

(17)  <F  f,f,>  >  - 1 -  e-(l-j)2/2a2  e-(l+J)2/2(2E2+a2)^ 

o  i  J  — ? — t— 

tt  a/2i  +CT 

2 

Figures  8  and  9  are  graphs  of  Eq(E)  for  various  truncated 
approximations.  The  two  left-arguments  of  the  function  CLERR  set 
the  number  of  neighbors  retained  in  the  particle  and  cell  trunca¬ 
tions,  respectively.  The  particle  mass  is  one  in  all  cases,  and 
a  is  set  equal  to  the  value  that  conserves  mass  in  the  truncated 
particle  approximation.  The  appearance  of  phenomenal  accuracy 
as  I  +  0  is  to  be  disregarded  since  it  refers  to  approximate 


I*-. 
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calculation  of  density  within  the  particle  at  X=0.  Figure  8 
approximates  the  untruncated  cell  over  the  range  0  <  z  <_  1 ,  and 
indicates  that  cell  error  is  limited  by  the  error  in  the  particle 
approximation.  Figure  9  shows  the  effect  of  truncation  of  the 
cell  approximation. 

The  final  figure  shows  1-p^  (with  m=l)  in  various  approxima¬ 
tions.  The  indication  is  that  mass  conservation  provides  a 
reasonable  estimate  of  a  value  of  E  for  which  truncation  error 
is  relatively  small  and  the  density  is  not  the  density  within  a 
particle.  In  the  case  where  the  number  of  neighbors  retained 
is  the  same  for  particles  and  cells,  the  estimate  is 

(IS)  l±  =  a1. 
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Difference  Analogs  of  Hamiltonian  Systems 

Frederic  Bisshopp 

Division  of  Applied  Mathematics 
Brown  University 
Providence,  Rhode  Island  02912 

Numerical  integration  of  Hamiltonian  systems  brings  forth 
the  question  of  whether  or  not  the  difference  analogs  have  an 
energy  invariant  that  corresponds  to  the  energy  integral  of 
the  differential  equations.  Several  difference  analogs  of 
the  harmonic  oscillator  will  be  examined  in  that  light,  and 
then  generalizations  will  be  given. 

The  differential  equations, 

•  • 

x  =  p,  p  =  -x,  (1) 

have  the  energy  integral, 

E  -  J(p2  +  x2) .  (2) 

Action-angle  variables  for  the  harmonic  oscillator  are  J=2irE 
and  0 ,  with 

x  =  /TE  cos(0-9o),  p  =  -  /^E  sin(0-0o)  (3) 

and 

0  =  1,  J  =  0.  (4) 

Even  the  simplest  difference  analog, 

ek  =  kh,  jk  =  jo,  (5) 


gives  exact  values  of  the  solution.  The  problem  of  interest 
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here  has  to  do  with  difference  analogs  of  the  primitive  equations 
(1)  • 

Explicit  difference  analog: 

xk+l  =  xk+hPk'  Pk+1  =  Pk‘hxk  (6) 

defines  the  iteration. 


where  the  reversed  arrow  indicates  new  values  that  are  assigned 
when  t  +-  t+h.  The  eigenvalues  are  1  ±  ih  and  their  magnitude 
is  /l+h2.  Thus  the  explicit  analog  is  unstable';  its  phase- 
portraits  spiral  outward. 

Implicit  difference  analog; 

xk+l-hpk+l  =  xk'  pk+l+hxk+l  =  pk  (8) 

defines  the  iteration, 


The  eigenvalues  are  (l±ih)/(l+h  )  and  their  magnitude  is 
l//l+h2.  The  implicit  analog  is  stable,  but  unsatisfactory 
since  its  phase-portraits  spiral  inward. 

Central  difference  analog: 

xk+l  "  1  hPk+l  =  xk  +  I  hPk 
pk+l  +  I  hxk+l  =  pk  "  I  hxk 


(10) 


defines  the  iteration 


l+jh 


(11) 


T  . 


The  matrix  that  transforms  (x,p)  is  a  rotation,  so 
12  2 

E  =  2^Pk+xk^  '''S  an  ener9y  invariant  of  the  central  difference 
analog.  The  only  source  of  truncation  error  is  in  the  phase 
of  the  solution. 


=  /2E  cos(kx-0o),  p^  =  -/2E  sin(kT-0o), 


(12) 


where 

tan  t=h/(l-  jh2).  (13) 

For  the  exact  solution,  x  would  be  h,  as  in  equations  (3,5); 
from  equation  (13) 

x  =  h  -  h3  +  °  (h5)  ,  | h |  <  2.  (14) 

The  attractiveness  of  the  central  difference  analog  is  diminished 
because  it  is  an  implicit  scheme. 

Leap-frog  method: 


xk+l  -  xk-l+2hpk'  pk+l  “  pk-l-2hxk 


(15) 


is  a  three-level,  explicit,  central  difference  scheme.  It  has 
solutions  of  the  form 


(16) 


where 


x2-l 


-2hX 


=  0  . 


x2-l 


2  2  /  7 

The  eigenvalues  are  solutions  of  X  =  l-2h  ±2ihvl-h  and  their 
magnitude  is  1  for  |h)  <  1.  There  is  a  well-behaved  mode, 

xk  =  /IE  cos(kT-0Q),  pk  =  -/2e  — j"~ ■  sin(kx-0o),  (18) 

and  an  ill-behaved  mode, 

x.  =  A(-l)kcos(ki-e  )  ,  p  =  A  -tJ—  (-1) ksin  (k T—  0  )  ,  (19) 

o  a  sin  t  o 


where 


tan  2t  =  2h  Jl-h2/ (l-2h2) . 


The  phase  error  follows  from 


t  =  h  +  \  h3  +  0  (h5)  ,  |  h  |  <  \/ft  , 


and,  if  the  ill-behaved  mode  is  suppressed,  the  energy  invariant 
of  the  leap-frog  .method  is  a  close  approximation  of  the  energy 
integral  since  h/sin  t  =  1+0 (h^). 


Flip-flop  integration: 

The  explicit,  two-level  scheme, 

xk+l  =  xk+hPk'  Pk+l+hxk+l  =  Pk 


defines  the  iteration. 


23 
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is 


The  eigenvalues  are  (2-h^±ih74-h  )/2,  and  their  magnitude 
1  for  |h|  £  2.  The  matrix  that  transforms  (x,p)T  clearly  is 
not  a  rotation,  but  an  energy  invariant  can  be  found  by 
multiplying  the  second  of  equations  (22)  by  (p^+^+p^) /2 . 

The  result. 


E  =  f<Pk+hpk*k+x2) 


i(U  +  !><Pk«k>2  +  (1  -  5)IPk-*k)2> 


(24) 


indicates  that  the  solution  lies  on  an  ellipse  at  the  angle 
0Q-kT,  where 

tan  t  =  hjl-h2/  (2-h2)  (25) 


The  phase  error  follows  from 

x  =  h  +  h3  +  0 (h5) ,  | h |  <  /2 


(26) 


By  comparison,  the  flip-flop  integration  is  best  at 

representing  the  phase  of  the  oscillation,  worst  at  representing 

the  energy  integral.  Nevertheless,  it  has  an  energy  invariant, 

2 

and  for  the  action-angle  variables  the  relative  errors  are  h  /4 
2 

and  h  /24,  respectively,  when  (xQ,po)  =  (xQ,0)  or  (0,pQ) .  By 
the  more  familiar  reckoning,  the  relative  error  in  the  waveform 
is  of  0(h). 

The  principal  advantages  of  the  flip-flop  integration  are 


that  it  is  two-level  and  explicit,  and  it  can  easily  be  generalized 
for  some  Hamiltonian  systems.  Consider  first 
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X  =  p,  p  =  f  (x) 
and  the  flip-flop  integration 


(27) 


xk+l  "  xk+hpk'  pk+l  pk  +  hf(xk+l}' 


(28) 


With  one  degree  of  freedom,  f(x)  is  derivable  from  the  potential 
x 


V  (x) 


f (z) dz. 


(29) 


and,  by  the  trapezoid  rule. 


V(*k+l,“V(V  =  ■  |<f(xk+i)+f(xk))  (X 

+  ij  *"(«)Uk+1-Xk)3 

where  E,  lies  in  (x^x^^).  We  define 
Ek  =  7  pk  "  7  hpkf(xk)  +  v(xk)f 


k+1 


-x 


k 


) 


(30) 


(31) 


and  it  follows  from  equations  (28,30)  that 

Ek+1-Ek  -  II  (32) 

It  cannot  be  seen  by  this  argument  whether  or  not  the  flip-flop 

integration  has  an  energy  invariant,  but  there  is  a  near 

invariant  that  confines  phase  portraits  of  solutions  of  the 

difference  equations  to  a  neighborhood  of  the  energy  integral, 

with  relative  errors  of  0(h)  for  times  as  long  as  0(l/h).  A 

2 

possible  drift  of  0 (h  )  for  times  of  0(1)  is  more  or  less 

2 

consistent  with  0 (h  )  errors  in  action-angle  variables,  and 


again  there  are  relative  errors  of  0 (h)  in  the  waveform  of  an 
oscillation.  It  can  be  expected  that  equation  (32)  overestimates 
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the  drift  per  cycle  of  a  nonlinear  oscillation. 

Now  let  us  consider  Hamiltonian  systems  defined  by 
H (x, p)  =  T(p)+V(x),  where  x  and  p  are  vectors.  Let 


-  I P=P  '  ^k  ~  —  I x=x 


(33) 


From  the  Taylor  series  of  T  and  V,  expanded  about  (xk,pk)  and 


*xk+l,pk+l^'  ifc  follows  that 


H(xk+l'pk+l)_H(xk'pk)  "  2(^k+gk+l)  (pk+l"pk) 


(34) 


2(fk+fk+l)  (xk+l~xk)+0 ( lxk+l_xkl  ''pk+l_pk'3) 


where  |£|2  =  £T£.  It  follows  easily  that  there  is  a  near 
invariant  (to  0(h3))  for  either  of  the  flip-flop  integrations. 


xk+l  =  xk+hgk  '  pk+l  =  pk+hfk+l' 

Ek  =  T(pk)+V(xk)  -  \  hgkfk, 
or 

pk+l  =  pk+hfk'  xk+l  =  xk+hgk+l  ' 

Ek  5  T(pk)+V(xk)  +  \  hgkfk. 

The  argument  fails  for  the  general  H(x,p). 

As  an  example,  the  cubic  nonlinear  oscillator, 

x  +  x^  =  0, 


(35) 


(36) 


(37) 


has  been  simulated  by  flip-flop  integration, 
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[1]  x  ■*-  x  +  hp 

[2]  p  «-  p  -  hx3  (38) 

[3]  -  1 

The  results,  as  illustrated  in  the  accompanying  figure,  suggest 
that  although  the  near  invariant  is  not  strictly  conserved 
there  may  be  no  systematic  drift  after  all.  The  phase  portrait 
of  approximately  twenty  cycles. of  the  oscillation  appears  to  be 
confined  to  a  band  (probable  width  of  0 (h  ) )  centered  about  the 

near  invariant.  The  value  of  h  in  the  computation  is  0.5  - 

for  most,  but  not  all  values  of  h  greater  than  a  number  near 
0.75,  the  algorithm  is  unstable. 


Figure  caption: 

Phase  portrait,  x  vs.  x  for  0  £  t  <^125  in  steps  of  0.5. 
Initial  data:  x(0)  =  1,  x(0)  =  0. 
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Appendix  C:  Plow  in  a  channel  by  biased  differences 

The  flow  that  is  simulated  here  occupies  a  rectangle  with 
nodal  coordinates 

Y  =  MB,  MM,  ...,  2,1  (bottom  to  top) 

X  =  1,  2,  NN,  NB  (left  to  right) 

The  sets  of  indices  that  are  used  to  define  central  nodes  and 
their  neighbors  are 


N1 

= 

1, 

...,  NN-1 

NO 

— 

2, 

. . .,  NN 

N3 

= 

3, 

. . . ,  NB 

m2 

= 

1, 

...,  MM-1 

MO 

= 

2, 

. . . ,  MM 

M4 

= 

3, 

. . . ,  MB 

Thus  the  correspondence  between  index-pairs  and  nodes  is 


Central  — 1  — 

(MO, NO) 

Western  - - 

(M0,N1) 

Northern  - - - 

(M2, NO) 

Eastern  - - ■ 

(M0,N3) 

Southern  — - - 

(M4,N0) 

The  input  variable,  u,  and  output,  z,  are  5  hy  MB  by  NB  arrays 
of  values  of  u-^,  u?,  f^,  f0  andirat  times  t  and  t+k...  The  parameter 
Nu  is  actually  vk/h2  since  length  and  time  are  measured  in  units  of 
h  and  k  to  eliminate  unnecessary  multiplications. 

The  first  part  of  the  program  carries  out  the  forward  integra¬ 
tion  of  X  =  u  as  described  in  Section  3  of  this  report.  The 
second  and  third  parts  are  decoupled  from  one  another  in  this 
early  version  of  the  method  because  only  the  hydrostatic 
component  of  the  pressure  boundary  conditions  has  been  included. 
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Thus  the  boundary  conditions  are: 

Western  boundary  (inflow) 

*c  =  %  '  hfl(^c) 
u  =  (constant,  0) 

— C 

Eastern  boundary  (outflow) 

tt  =  constant 
c 

u0  =  (ulM.0) 

Northern  boundary  (rigid) 

*c  =  "s  +  hf2^> 

“c  =  (0,0) 


(hydrostatic ) 
(slug  flow) 


(fully  developed) 


(hydrostatic ) 
(no  slip) 


Southern  boundary  (center  line) 

”c  -  "n  -  hf2(2c) 

uc  -  (uln,0) 


(hydrostatic ) 
( symmetric ) 


The  last  part  of  the  program  was  added  to  correct  for  global 
truncation  error.  At  each  station  along  the  x-^-axis  u^  is 
renormalized  to  preserve  the  constant  total  mass  flux. 

The  initial  condition  was  slug  flow  everywhere  but  at  the  rigid 
boundary  of  an  11  by  31  node  grid.  The  initial  value  of  Uj 
and  the  inflow  velocity  were  set  at  1A/2  to  insure  that  no 
bias  lines  would  fall  outside  local  skeletons,  and  the  body 

p 

force  was  set  to  zero.  Nu  (=  vk/h  )  was  set  at  0.2,  and  after 
15  time  steps  (CPU  time,  12  sec.)  the  flow  was  almost  at  a  steady 
state.  Figure  1  shows  the  profiles  of  u^  at  the  horizontal 
stations  1,  3,  5,  7  and  9-  The  numbers  above  and  below  each  graph 
are  the  maximum  and  minimum  values  of  u^.  In  figure  2  the  same 
profiles  are  shown  after  10  more  time  steps.  In  figures  3  ~  5 
the  value  u?  =  0  at  the  rigid  boundary  is  not  included  and  the 
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scale  is  expanded  to  show  detailed  structure  of  the  velocity 
profiles.  The  progression  of  the  point  of  maximum  velocity 
from  edge  to  center  as  one  moves  downstream  is  consistent  with 
results  of  other  numerical  simulations  and  with  results  of 
boundary  layer  theory. 
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FIGURE  2; 


VELOCITY  PROFILES 

STATIONS  1,3, 5, 7, 9  /1PTEP  2  5  TlttE  STEPS 
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FIGURE  3:  STATIONS  1  TO  8  AFTER  2  5  TIME  STEPS 
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Appendix  D:  Thermal  instability 

The  problem  that  was  used  to  check  the  biased  difference 
algorithm  for  thermal  convection  was  the  B^nard  problem.  Thermal 
instability  of  a  horizontal  layer  of  a  fluid,  heated  from  below, 
is  governed  by  the  equations 

X  =  u 

u  =  -  Ejt  +  agz(T  -  <T>)  +  vAU 

V/u  =  0 

T  =  kAT 

where 

if  =  -2-  +  gz  , 

Po 

a  is  the  coefficient  of  thermal  expansion,  g  is  the  acceleration 
of  gravity,  v  is  the  kinematic  viscosity,  k  is  the  thermal 
diffusivity,  and  <T>  is  the  mean  temperature.  The  boundary 
conditions  that  were  used  are 

0  all  boundaries 

Tq  and  T-^  lower  and  upper  boundaries 

0  sidewalls 

Two-  and  three-dimensional  cases  were  implemented.  For 
simplicity  we  will  discuss  the  two-dimensional  case;  the  three- 
dimensional  algorithm  is  Identical  in  form.  Let  the  subscripts 
1,  2  denote  the  components  of  x  =  (y,z),  and  let  the  subscripts 
c,  n,  s,  u,  d  denote  central,  northern,  southern,  upward,  and 
downward  nodes  in  a  rectangular  grid.  In  the  APL  function  TW1C 
the  4  by  M  by  N  dimensional  array  U  contains  nodal  values  of 
u,  T  and  it  .  North  is  to  the  right  and  up  is  up  in  the  arrays; 
the  correspondence  between  index  pairs  and  directions  is 


u  = 

T  = 
3T 

3n  = 


•  ■'•i.jm.  m. 
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U  « — U[;MO;NO] 

v 

u  * — »  U  [ ;  MO ;  N1  ] 

S 

Uu  «« — U[;M2;N0] 
Un  < — U[;M0;N3] 
Ud  *  U [ ; M4 ; NO ] 


N1  =  NO-1 
M2  =  MO-1 
N3  =  Nl+2 
M4  =  M2 +2 


The  first  part  of  the  algorithm  takes  input  values  u,  T,  irat 
time  t  (say);  the  parameters  P  =  v,  k,  ag  are  not  used  until 
later.  The  forward  integration  of  X  =  u  with  time  increment 
k  gives 

xx  =  +  ku( xc ) 


where  xc  is  the  position  at  time  t  of  the  bias  line  that  passes 

through  xc  at  time  t+k,  and  u  still  refers  to  time  t.  With 

linear  interpolation  for  u(x )  the  displacement 

c 


is  given  by  the  linear  equations 
(I  +  kM)6.  =  -  ku 

V 

where  the  space  increment  is  h  and 

M11  =  K«uln"  ulc^ulc  >0)  +  <ulc"  ulsXulc  <0)) 

M12  =  h^ulu"  ulc^u2c  >0^  +  (ulc"  uld^u2c  <0^ 

M21  =  h^u2n"  u2c^ulc  >0^  +  (u2c_  u2s^ulc  <0^ 

M22  =  H^u2u"  u2c^u2c  >0)  +  (u2c“  u2d^u2c  <  °^  * 

Note, the  parenthetic  expressions,  (uic<0)  and  (uiC >  which 
have  logical  values  0  or  1,  are  used  to  compute  arrays  of 
coefficients  of  M  in  the  APL  function. 


Given  arrays  of  S_  at  the  central  nodes,  the  initial  data  for 
velocity  and  temperature  is 

u  =  u  +  MS 
— c  — c  — 

T  =  T  +  m*<S 
c  c - 

where  the  vector  m  has  components 

">1  *  K<<Tn  '  Tc><ule  >0>  +  <Tc  -  Ts»ulo  <°>> 
m2  =  K«\,  -  T0)(U20  >0>  +  ‘  Trl)(U5^  <0D- 


‘c  d ' '  2c 


At  this  point  the  approximation 

(Z‘u)c  =  2h^Uln  “  uls  +  U2u  “  U2d ^ 

is  evaluated  for  use  in  the  pressure  equation. 

The  backward  integration  of  the  temperature  equation  defines 
the  iteration  that  initiates  the  second  part  of  the  algorithm 

Tc  =  (fc  +  £|<Tn  -  Ts  +  Tu  +  Ta))/(1  +  *£)  . 


Substitution  of  =  0  in  the  divergence  of  the  backward  time 
Integration  of  the  momentum  equation  gives  the  pressure  equation 
also  iterated 

%  =  *(7Tn  +  "s  +  uu  +  ff-d  +  Td>  -  T-(r*)c)  • 


The  momentum  equations, 


ulc  ~  ^ulc  +  2h(lrs“  ^  +  ^uln+  Uls+  Ulu+  Uld^^1+  ^2^ 

u2c  =  [Z2c  +  Ih^d-  +  g^(Tc-<T» 

+  r|(u2n+  u2s  +  u2u+  u2d^/^1  +  TT^ 
n  n 


complete  the  algorithm  at  internal  nodes. 

The  temperature  is  fixed  (unchanged)  at  the  top  and  bottom 

boundaries,  T  =  T  at  the  south  boundary  and  T  =  T_  at 
c  n  c  s 

the  north.  Both  components  of  u  are  zero  on  all  boundaries, 
and  the  pressure  boundary  conditions  are 

Top: 

"c  =  7rd+  5-aSh(TG+  Td“  2<T>)  +  ^-u2d+  Tj^(u2d~  u?d) 

Bottom: 

77 c  =  17 u“  5txS^(Tc+  Tu-  2<T>)  -  u2u  -  t^(u2u-  u2u) 

North: 

,  2v  ,  h  /  —  \ 

"c  =  *8  +  "PT^lS  +  ^Uls  "  Uls) 

South: 

77  c  =  Tn  ‘  ifSn  '  W^ln  "  uln^ 

The  small  terms  of  -^(u  -  u)  were  not  included  in  the  versions 
of  the  algorithm  that  were  used  here. 

Finally,  the  last  part  of  the  algorithm  renormalizes  horizontal 
and  vertical  components  of  velocity  to  account  for  V-u  =  0 
globally,  and  sets  the  undetermined  constant  in  the  pressure 
field  at  a  convenient  value. 

The  onset  of  thermal  instability  in  a  fluid  layer  of  depth 
d  occurs  at  sufficiently  high  values  of  the  Rayleigh  number 

R  _  ctgATd^ 

KV 

where  AT  is  temperature  at  bottom  minus  temperature  at  top. 

For  an  infinite  layer  bounded  above  and  below  by  rigid  boundaries 
the  critical  Rayleigh  number  is  1708,  and  for  R  greater  than 
that  convection  appears.  The  fluid  motion  is  in  vortices 


with  alternating  clockwise  and  counterclockwise  circulation 
(rolls),  and  the  width  of  a  vortex  is  very  slightly  greater 
than  its  depth,  d. 

In  the  simulations  we  have  run,  time  and  distance  are 
measured  in  units  of  k  and  h,  i.e.  At  =  Ax  =  1.  Thus  an 
M  by  N  array  of  nodes  has  a  width  (N-l)  and  a  depth  (M-l), 
and  the  Rayleigh  number  is 

R  _  ctgAT(M  -  l)3 

KV 

For  the  two-dimensional  cases  the  grid  is  11  by  30,  k  =  v  -  AT  =  1, 
and  ag  is  set  at  2  or  4  to  give  Rayleigh  numbers,  2000  or  4000. 
Runs  were  started  with  a  uniform  temperature  gradient  and 
small,  randomly  chosen  velocity  perturbation.  Features  of  the 
resulting  steady  states  that  developed  are  shown  in  Figures  1-6. 

In  figures  1  -  3,  R  is  2000  and  there  are  two  weak  primary 
cells  and  a  very  weak  secondary  cell.  This  is  in  accord  with 
stability  theory  which  predicts  that  neither  two  nor  three  cells 
can  fill  a  container  with  aspect  ratio  3:1  at  slightly  super¬ 
critical  Rayleigh  numbers.  As  far  as  we  know,  there  is  no 
theoretical  work  that  has  predicted  that  the  resolution  of  that 
would  be  a  secondary  cell  that  is  driven  by  the  primary  motion 
rather  than  by  buoyancy. 

In  figures  4  -  6,  R  is  4000  and  no  secondary  cell  is  needed  to 
fill  the  container.  This,  incidentally,  is  about  the  limit  of 
what  can  be  simulated  with  At  =  Ax.  When  a  component  of  uk/h 
exceeds  one,  the  bias  line  falls  outside  the  local  skeleton, 
interpolations  become  extrapolations,  and  the  method  begins  to 
be  unreliable  and  possibly  unstable. 

Finally,  a  three-dimensional  cell  was  simulated  at  Rayleigh 
number  4000  in  an  11  by  11  by  11  node  container.  The  alignment 
of  the  circulation  with  its  axis  on  a  diagonal  in  plan  view 
allows  a  cell  that  has  a  width  that  is  slightly  greater  than 
its  height.  Figure  7  is  a  plan  view  of  the  contours  of  w  at  z  =  6 
accompanied  by  elevation  views  of  the  direction  fields  of  u 
and  w  at  y  =  6  (below)  and  v  and  w  at  x  =  6  (to  the  left  ). 


Altogether,  biased  differences  appear  to  be  a  reliable  and 
effieibnt.  way  to  simulate  thermal  convection.  CPU  time  on 
the  IBM  370/158  at  Brown  University  is  5  seconds  per  time  step 
of  the  three-dimensional  simulation  on  1331  nodes. 
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FIGURE  2:  R  = 2000 

VERTICAL  VELOCITY  AT  MIDHEIGHT 
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FIGURE  3:  R= 2000 

TEMPERATURE  AT  MID  HEIGHT 
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FIGURE  5 :  R  = 4000 

VERTICAL  VELOCITY  AT  MIDHEIGHT 
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FIGURE  6:  i?=4000 

TEMPERATURE  AT  MID  HEIGHT 
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Appendix  F:  Interactive  generation  of  contour  maps 

The  APL  function  CONTOUR  will  be  described  here.  It  takes 
as  Input  a  rectangular  array  F  of  real  numbers  that  can  be 
regarded  as  heights  above  a  rectangular  grid  of  nodes  In  the 
x-y  plane.  The  second  input,  H,  is  one  or  several  values  for 
which  one  wants  the  level  lines.  The  output,  X,  is  a  sequence 
of  coordinates  X,Y,X,Y, . . . ,-l,-l,X,Y, . . .  of  points  on  the  level 
lines.  A  pair  of  entries  -1,-1  within  X  signals  the  start  of  a 
new  contour  line. 

As  always,  there  is  the  decision  to  be  made  whether  the  array 
F  represents  nodal  values  of  a  function  F(X,Y)  or  values  of  H 
at  positions  in  the  array.  CONTOUR  does  neither  case:  If  F  is 
nodal  values  of  F(X,Y)  one  must  first  take  the  transpose  of  F, 
and  if  F  is  values  of  H  above  an  array  one  must  first  reverse 
the  order  of  the  rows  of  F.  Thus 

H  CONTOUR  $  F  (function) 

or 

H  CONTOUR  0  F  (array) 

gives  results  in  standard  position  with  X  Increasing  to  the  right 
Y  increasing  upward. 

Contour  map  algorithms  require  a  check  for  F  =  H  at  a  node  of 
the  array  because  the  contour  can  have  one  or  more  branches  there 
(e.g.  at  a  saddle  point).  The  interactive  aspect  of  the  function 
operates  when  such  a  point  is  found,  and  the  user  is  required  to 
enter  a  perturbed  value  of  H.  Then,  at  a  saddle  for  example. 
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The  function  CHECK,  called  by  CONTOUR,  first  checks  F  =  H 
at  corners  (nodes  of  the  array);  then  introduces  centers  of 
the  rectangular  patches  where  F  is  assigned  the  average  of  the 
values  at  corners;  then  checks  F  =  H  at  centers.  Given  F  /  H 
at  any  corner  or  center  of  the  now  triangulated  array,  no 
intersection  or  branching  of  contour  lines  is  possible,  and 
the  only  decision  left  is  whether  a  contour  is  closed  and  lies 
within  the  array  or  open  and  terminates  at  an  edge. 

The  function  NODES  finds  nodes  on  the  contour  by  linear 
interpolation  of  values  of  F  at  corners  and  centers.  The 
output  of  NODES  has  four  rows  and  as  many  columns  as  the 
number  of  nodes  on  the  contour.  The  first  two  rows  are  encoded 
addresses  of  line  segments  that  have  a  node,  and  the  last  two 
rows  are  coordinates  of  the  nodes.  Finally,  the  function  SORT 
uses  the  encoded  addresses  to  arrange  the  nodes  in  order,  checks 
for  open  or  closed  contours,  and  returns  the  coordinates. 

The  logical  structure  of  this  algorithm  is  much  .simpler  than 
others  we  have  examined,  and,  even  in  APL,  it  is  considerably 
faster  than  its  FORTRAN  competitors. 
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□Ci? ' CONTOUR ' 

A>/?  CONTOUR  F  \N \S  ’,E ’,W iUIO 
A  X~X, Y,X, . . .X, X,Y 
ft  ON  CONTOURS  OF  H=<HF(X.Y)  OR  H=eFiI  -,J3 
A  CALLS  VX*- CHECK  HV 

ft  INDICES  -  NORTH, SOUTH,  EAST,  WEST  ON  e F 

nio*-o 

N*-  1+S+  »(pF)[  Oj-1 
E«-l+J/-<-t(pF)[l]-l 
X+~2  *  CHECK  H 


UCR '  CHECK ' 

X-*- CHECK  H\Z\ZC 

A  EXTERNAL  VARS  F ,N,S,  E,  W,QI0  =  Q 
A  CALLS  VX+SORT  JV  A//D  VY+NODESV 
a  CHECK  AND  GO  -  ANY  NUMBER  OF  CONTOURS 
-►(  1  =p  ,U)/BCN 
JM  CHECK  Hi  0  2),  CHECK  HH 
-*•0 

BCNiZ+F-H 

A  CONTOUR  THROUGH  A  CORNER 
+(*/ */UCT<\Z)/OKl 

O-'Fs*  ,(1H),  '  AF  A  CORNER  -  TRY  //=' 

//•♦•aB 

■*BGN 

A  CONTOUR  THROUGH  A  CENTER 

OKI :ZC+0 .25* ,+/Z[S iW3 ,Z[S iEl ,ZL NiEl, Cl. 5] 

-+(a/L)CT<IZC)/OK2 

CK'F  =  ’ ,(▼//),  '  A?  A  CENTER  -  27?/  //=' 

//«-aB 

■*BGN 

A  WO  CONTOUR 

OK2i-+(v/v /b*\+/xZlS-,W2,ZiS-,E2  ,ZlN-,E)  ,[1. 5]  ZCffjJ/]  )/CF3 
EK'tfC  CONTOUR  FOR  H  =  ',{lH),'  -  TRY  H=' 

H «  aB 
■+BGN 

OK3  :  X+-S0RT ANODES 


•i  J* 
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□Cif '  NODES » 

Y<rNODES-,K-,X-,V 

A  EXTERNAL  VARS  N ,S ,E ,U ,Z ,ZC ,010=0 
A  ADDRESSES  AND  POSITIONS  OF  NODES 
4  o  po 

A  N-S  NODES 

NS  :  ■+(  0=pX^(  ,  )*xz[S;]  )/i(  l  +  pE)*p  N) /EW 

^(6xx).^^((pW).i+pE)ta: 

Y+X,<*X--3,Q,(.  y+(  ,Z[//;])C^]-^(,ZC5i])[X]),C0.  5]  0 
A  E-W  NODES 

EU '.-*•(  0  -pK*-  (  ,  (  *ZC  iE2  )*x£[  ;  f/  J  )  /  i(  p  E)*l  +  pN )  / SU 
X+(6*X),X+6a(l+pN)  ,pE)7K 

X+X ,$X-Q,~3t0,l0 .5]  V*(  .ZliEl)lia-V+(  ,Zl  iWDlKl 

A  SJ/-C  JVOD&i’ 

0=pK*-UxZC)*,*ZlSiWl)/  \  (.p  E)*pN)  /SE 
X+(6xX),X+$UpN).pE)rK 

X-*-X ,  QX-~2 ,  ~2  ,  V,  [  0 . 5  ]  V+-V*2xZCLK]-V*-(  ,Z[SiN])LKj 
A  SE-C  NODES 

SE  :  -*■(  0  =  pA'«-((  xZC)*,  xZlS-.El )/  i(p£)xpif)//if£ 
^(6x/),^((p/)p  o  l)+X*-t)((pN),pE)TX 
X-X.<HX-~2.2,  V,[0.S]-V+V*2xZC[K]-V+(.Z[SiE])[X] 
a  NODES 

NE:-(Q=pK<-((  xZC)*.xZlUiEl)/ \(pE)xpN)/NW 
X+(6xX),X+-l+<*((pN)  ,pE)TK 

X+X,<HX+~2,~2,V,lO.  5]  V+Vi2xZCUQ-V+(  ,Z[ //;£•])[ K] 
a  AW- <7  iWPFS 

NW  0  =pa'-*-  ((  xZ  C) * ,  xZt//;V]  )/i(p£)*p//)/0 

X+(&xX)  ,X+((pX)p  1  0  )  +X+-§(,  ( p  N ) ,  p  E  )tK 

X*-X  ,<HX+~2, 2,  K,  to.  5]-^7*2xZCC^3-^(  ,ZCW;f/]  )CA:] 


□  CR' SORT' 

X+SORT  Y \N \L\M\K  \ I 
A  EXTERNAL  VAR  DIO=0 
X+iO 

BGN:N/+-iL*-(pY)[  0] 

a  START  WITH  ZERO  AND  ONE  OR  TWO  NEIGHBORS 
N*-~li>0,  (  V/  (+/  \IlLp  Oil-I+Yl ;  0, 1 J  )  •.  =  2  3)/// 

A  CHECK  FOR  ENDPOINT  OR  CLOSED  CON  TO  \H 
B:-+{  0  =pA'+ (£*;/[  1]) /A' «-(v/(*/  |  J [ £,p //[  Oj  ;  = 

■+{K*~lU1<rK,N)/B 
■*END 
C :  N+-4>N 

a  FIND  OTHER  ENDPOINT 

DiN+(.K-(K*Nlll  )/A^(v/(  +  /  |  JCLpA'C  0]  ;]-!)•.  -  2 
-*■(  0*pA)  ID 

END:X<-X,(  ,y[A;3,2]),2p~l 
-*-(LspiV)/0 

-+BGN 


UPPLEMENTARY 


INFORMATION 


p&> 


SECURITY  CLASSIFICATION  OF  THIS  PACE  (Wh» n  D.l.  F.nertd) 


r\Kjj  /u\j7~o*  o** 


VV. 


1  REPORT  DQCUMEHTATlOH  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  []  / ]  fj 

2)  30VT  ACCESSION  NO. 

J.  RECIPIENT'S  CATALOG  NUMBER 

A  title  f«r iB  SuhUII.J 

Final  report  on  manuscripts 
submitted  or  published  under  the  subject 
Grant. 

_ 

S.  TYPE  OF  REPORT  A  PERIOD  COVEREO 

Final  Report 

15  June  76  -  15  Sept.  J) 

6.  PERFORMING  ORG  REPORT  NUMBER 

1  AUTHORf.)  f 

““ _ 0 

6.  CONTRACT  OR  GRANT  NUMBERfaJ 

DAAG29-76-G-0178^ 

s 

».  PERFORMING  ORGANIZATION  NAME  AND  AODRESS  _ ^ 

Division  of  Applied  Mathematics 

Brown  University 

Providence,  R.I.  02912 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  6  WORK  UNIT  NUMBERS 

Project  No.  1354m 

11.  CONTROLLING  OFFICE  NAME  AnO  ADDRESS 

'2  REPOET  DATE 

U.  o.  Army  Research  Office 
-  tm-  Bex  -L-£r’P  1 1  - 


12  July  1930 

’3  NUMBER  OF  PAGES 


Research  Triangle  Park,  27709 

Ti  MONI  TOPING  AGENCY  NAME  *  AOORESSflf  dlltmrfU  tram  Crnlrolllnt  OWct)  IS 


4 

SECURITY  CLASS,  fof  (hit  rtport) 


16 


Ur. c lacs i  fied 

is*  DECL  A«TiF|CATIOn/OOWNGRADIN 
schedule 


DISTRIBUTION  STATEMENT  fof  IM<  Rtport) 


"~1 


Approved  for  public  release;  distribution  unlimited. 


*7  DISTRIBUTION  STATEMENT  (of  the  ebatrect  entered  In  Block  20,  It  different  from  Report) 


1®  SUPPLEMENT  ARV  NOTES 


19 


The  view, 
author1,  ) 
positi  ;r. , 


pini  cr.s , 
arid  should 
po 1 i ey  ,  r 


ar.  i/or  findings  c 
n  t  be  construed 
Ivcision,  unless 


1  'n ta i ned  in  this  report  are  those  of  the 
as  an  official  department  of  the  Army 
so  lerignated  by  other  documentation. 


KCv  WOROS  (Continue  on  revere*  aide  ft  neceeemry  end  Identity  by  block  number) 


20  ABSTRACT  (Continue  on  reveree  elde  It  neceeemry  end  Identify  by  block  number) 


00 


roRM 

1  J AN  7| 


1473 


EDITION  OF  I  NOV  6S  IS  OBSOLETE 


\ 


J 


rclassi fj ed 


FINAL  REPORT 


I 

1.  ARO  PROPOSAL  NUMBER:  13548M 

2.  PERIOD  COVERED  BY  REPORT:  15  JUNE  1976  -  15  SEPT.  1979 

3.  TITLE  OF  PROPOSAL: 

FINITE  ELEMENT  METHODS  FOR  HEAT  TRANSFER  PROBLEMS 

4.  CONTRACT  OR  GRANT  NUMBER:  DAAG29-76-G-0178 

5.  NAME  OF  INSTITUTION:  BROWN  UNIVERSITY 

6.  AUTHOR  OF  REPORT:  FREDERIC  BISSHOPP 


1354M 


DR.  FREDERIC  BISSHOPP 

DR.  BRUCE  CASWELL 

BROWN  UNIVERSITY 

DIVISION  OF  APPLIED  MATHEMATICS 

PROVIDENCE,  R.  I.  02912 


LIST  OF  MANUSCRIPTS  SUBMITTED  OR  PU  ILISHED  UNDER  ARO 
SPONSORSHIP  DURING  THIS  PERIOD: 

by  FREDERIC  BISSHOPP 


1.  A  Lagrangian  Finite  Element  Method  for  Barotropic  Fluids, 
(long  version)  This  will  be  rewritten  and  submitted  for 
publication  when  the  numerical  work  is  completed. 

2.  A  Lagrangian  Finite  Element  Method  for  Barotropic  Fluids, 
(short  version)  This  was  presented  at  the  International 
Symposium  oh  Innovative  Numerical  Analysis  in  Applied 
Engineering  Science,  Versailles  (France) ,  23-27  May  1977. 

3.  Difference  Analogs' of ^Hamiltonian  Systems. 

This  will  be  submitted  for  publication  upon  completion 
of  an  application  of  the  method  that  gives  an  efficient 
algorithm  for  construction  of  particle  paths  in  two- 
dimensional  flows. 

4.  Biased  Differences  for  Advection/Dif fusion  Problems. 

(long  version)  Three  working  papers  on  this  subject  will 
be  rewritten  as  one  or  two  papers  to  be  submitted  for 
publication. 

5.  Biased  Differences  for  Advection/Dif fusion  Problems. 

(shoit  version)  This  was  presented  at  the  Gordon  Con¬ 
ference  on  Fluids  in  Permeable  Media,  Plymouth,  New 
Hampshire,  June  1978. 

6.  Computer  Graphics  Software. 

This  is  now  available  in  the  APL  software  library  at 
Brown  University. 


by  MARION  MICHAUD 

Numerical  Simulation  of  Reservoirs. 

Ph.D.  dissertation,  August  1979 

Accurate  Waterflood  Simulation  Using  Biased  Differencing 
and  Selective  Grid  Refinement. 

Latest  draft,  March  1980.  To  be  submitted  for  publication. 


7. 


Continued  from  previous  page. 


by  E.  W.  FLERI 


1.  A  Lagrangian  Finite  Element  Method. (probable  title) 
Ph.D.  dissertation  in  preparation 


7a.  MANUSCRIPTS  SUBMITTED  OR  PUBLISHED  DURING  THIS  PERIOD 
WITH  PARTIAL  SUPPORT  UNDER  ARO  SPONSORSHIP: 


by  W.  K.  LYONS 


1.  The  Single  Conservation  Law  of  Discontinuous  Media. 

Ph.D  dissertation.  May  1980 

by  T.  G.  McKEE,  Jr. 

1.  Instability  and  Bifurcation  for  Two-Dimensional  Regions  of 
Constant  Vorticity. (possible  title) 

Ph.D  disseration  in  preparation. 


8.  SCIENTIFIC  PERSONNEL  SUPPORTED  BY  THIS  PROJECT  DURING 
THIS  REPORTING  PERIOD: 


6/76  to  5/77 

Summer 

Academic  Year 

Bisshopp,  F. 

2  mo. 

9  mo.  (10%) 

Fieri,  E.  W. 

— 

9  mo . 

Michaud,  M. 

1  mo. 

6/77  to  5/78 

Bisshopp,  F. 

2  mo. 

9  mo.  (10%) 

Fieri,  E.  W. 

2  mo. 

9  mo. 

Michaud,  M. 

1.5  mo. 

6/78  to  8/79 

Bisshopp,  F. 

1  mo. 

9  mo. (10%) 

Fieri,  E.  W. 

3  mo. 

— 

Michaud,  M. 

3  mo. 

5  mo . 

McKee,  T.  G 

2  mo 

9  mo. 

Lyons,  W.  K. 

— 

5  mo. 

w 


8a.  DEGREES  GRANTED  DURING  THIS  REPORTING  PERIOD: 

MARION  CATHARINE  MICHAUD,  Ph.D.,  June  1980. 

WILLIAM  KIMBEL  LYONS,  Ph.D.,  June  1980 

8b.  DEGREES  EXPECTED  FOR  WORK  SUPPORTED  BY  ARO: 

E.  W.  FLERI ,  Ph.D.  (June  1981) 

T.  G.  McKEE, Jr.,  Ph.D.  (June  1981) 

9.  RESEARCH  FINDINGS 

A  two-volume  technical  report,  dated  2  April  1980  and 
submitted  to  the  U.  S.  Army  Research  Office,  contains  brief 
descriptions  of  our  research  on  numerical  fluid  mechanics 
in  its  first  seven  sections.  More  detailed  descriptions  of 
parts  of  that  research  are  provided  in  six  accompanying 
appendices. 


